# Need to install dolfyn, easy to do it with: pip3 install dolfyn
# Ivica, 3-May-2022
# Import core DOLfYN functions
import dolfyn as dlfn
# Import ADCP-specific API tools
from dolfyn.adp import api
from matplotlib import pyplot as plt
import matplotlib.dates as dt
import numpy as np

ds = dlfn.read('/home2/ivica/notebooks/ADCP/ADCP1000.000')

Reading file /home2/ivica/notebooks/ADCP/ADCP1000.000 ...

ds

<xarray.Dataset>
Dimensions:              (beam: 4, range: 40, time: 11302, dir: 4, x: 4, x*: 4,
                          earth: 3, inst: 3)
Coordinates:
  * beam                 (beam) int64 1 2 3 4
  * range                (range) float64 3.09 4.09 5.09 ... 40.09 41.09 42.09
  * time                 (time) datetime64[ns] 2021-08-16T05:00:00 ... 2021-1...
  * dir                  (dir) <U3 'E' 'N' 'U' 'err'
  * x                    (x) int64 1 2 3 4
  * x*                   (x*) int64 1 2 3 4
  * earth                (earth) <U1 'E' 'N' 'U'
  * inst                 (inst) <U1 'X' 'Y' 'Z'
Data variables: (12/19)
    number               (time) uint32 1 2 3 4 5 ... 11299 11300 11301 11302
    builtin_test_fail    (time) bool False False False ... False False False
    c_sound              (time) float32 1.553e+03 1.553e+03 ... 1.518e+03
    depth                (time) float32 40.0 40.0 40.0 40.0 ... 40.0 40.0 40.0
    pitch                (time) float32 -0.31 -0.31 -0.3 ... -1.26 -1.26 -1.26
    roll                 (time) float32 2.5 2.49 2.48 2.48 ... -3.14 -3.14 -3.14
    ...                   ...
    vel                  (dir, range, time) float32 -0.618 -0.247 ... nan nan
    amp                  (beam, range, time) uint8 84 85 85 85 ... 154 161 147
    corr                 (beam, range, time) uint8 68 67 67 66 ... 77 82 90 82
    prcnt_gd             (beam, range, time) uint8 28 12 11 21 14 ... 0 0 0 0 0
    beam2inst_orientmat  (x, x*) float64 1.462 -1.462 0.0 ... -1.034 -1.034
    orientmat            (earth, inst, time) float64 0.3861 0.3858 ... -0.9983
Attributes: (12/39)
    inst_make:                TRDI
    inst_model:               Workhorse
    inst_type:                ADCP
    rotate_vars:              ['vel']
    has_imu:                  0
    name:                     wh-adcp
    ...                       ...
    water_ref_cells:          [1, 5]
    fls_target_threshold:     50
    xmit_lag_m:               0.5
    fs:                       0.004545454545454545
    vel_gps_corrected:        0
    h_deploy:                 0.6

xarray.Dataset

• Dimensions:
  • beam: 4
  • range: 40
  • time: 11302
  • dir: 4
  • x: 4
  • x*: 4
  • earth: 3
  • inst: 3
• Coordinates: (8)
  • beam
    (beam)
    int64
    1 2 3 4

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

  • range
    (range)
    float64
    3.09 4.09 5.09 ... 41.09 42.09

    units :

    m

    array([ 3.09,  4.09,  5.09,  6.09,  7.09,  8.09,  9.09, 10.09, 11.09, 12.09,
           13.09, 14.09, 15.09, 16.09, 17.09, 18.09, 19.09, 20.09, 21.09, 22.09,
           23.09, 24.09, 25.09, 26.09, 27.09, 28.09, 29.09, 30.09, 31.09, 32.09,
           33.09, 34.09, 35.09, 36.09, 37.09, 38.09, 39.09, 40.09, 41.09, 42.09])

  • time
    (time)
    datetime64[ns]
    2021-08-16T05:00:00 ... 2021-12-...

    array(['2021-08-16T05:00:00.000000000', '2021-08-16T05:15:00.000000000',
           '2021-08-16T05:30:00.000000000', ..., '2021-12-11T21:45:00.000000000',
           '2021-12-11T22:00:00.000000000', '2021-12-11T22:15:00.000000000'],
          dtype='datetime64[ns]')

  • dir
    (dir)
    <U3
    'E' 'N' 'U' 'err'

    ref_frame :

    earth

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

  • x
    (x)
    int64
    1 2 3 4

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

  • x*
    (x*)
    int64
    1 2 3 4

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

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

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

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

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

• Data variables: (19)
  • number
    (time)
    uint32
    1 2 3 4 ... 11299 11300 11301 11302

    units :

    array([    1,     2,     3, ..., 11300, 11301, 11302], dtype=uint32)

  • builtin_test_fail
    (time)
    bool
    False False False ... False False

    units :

    array([False, False, False, ..., False, False, False])

  • c_sound
    (time)
    float32
    1.553e+03 1.553e+03 ... 1.518e+03

    units :

    m/s

    array([1553., 1553., 1553., ..., 1518., 1518., 1518.], dtype=float32)

  • depth
    (time)
    float32
    40.0 40.0 40.0 ... 40.0 40.0 40.0

    units :

    m

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

  • pitch
    (time)
    float32
    -0.31 -0.31 -0.3 ... -1.26 -1.26

    units :

    deg

    array([-0.31, -0.31, -0.3 , ..., -1.26, -1.26, -1.26], dtype=float32)

  • roll
    (time)
    float32
    2.5 2.49 2.48 ... -3.14 -3.14 -3.14

    units :

    deg

    array([ 2.5 ,  2.49,  2.48, ..., -3.14, -3.14, -3.14], dtype=float32)

  • heading
    (time)
    float32
    112.7 112.7 112.7 ... 142.4 142.4

    units :

    deg

    array([112.72, 112.7 , 112.7 , ..., 142.41, 142.41, 142.41], dtype=float32)

  • temp
    (time)
    float32
    31.76 31.76 31.76 ... 17.27 17.26

    units :

    C

    array([31.76, 31.76, 31.76, ..., 17.27, 17.27, 17.26], dtype=float32)

  • salinity
    (time)
    float32
    38.0 38.0 38.0 ... 38.0 38.0 38.0

    units :

    psu

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

  • min_preping_wait
    (time)
    float32
    0.07 0.07 0.07 ... 0.07 0.07 0.07

    units :

    s

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

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

    units :

    deg

    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 :

    deg

    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 :

    deg

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

  • vel
    (dir, range, time)
    float32
    -0.618 -0.247 -0.243 ... nan nan

    units :

    m/s

    array([[[-0.618, -0.247, -0.243, ...,  0.038,  0.083,  0.108],
            [-0.418,  0.073, -0.067, ...,  0.037,  0.062,  0.099],
            [   nan,    nan,    nan, ...,  0.034,  0.069,  0.105],
            ...,
            [   nan,    nan,    nan, ..., -0.511, -0.563, -0.603],
            [   nan,    nan,    nan, ..., -0.288, -0.49 , -0.226],
            [   nan,    nan,    nan, ..., -0.212, -0.298, -0.354]],
    
           [[ 0.064,  0.02 ,  0.013, ...,  0.075,  0.09 ,  0.103],
            [ 0.165,  0.066,  0.088, ...,  0.054,  0.09 ,  0.081],
            [   nan,    nan,    nan, ...,  0.032,  0.077,  0.056],
            ...,
            [   nan,    nan,    nan, ..., -0.499, -0.487, -0.706],
            [   nan,    nan,    nan, ..., -0.364, -0.403, -0.429],
            [   nan,    nan,    nan, ..., -0.295, -0.321, -0.324]],
    
           [[-0.086, -0.04 , -0.04 , ..., -0.001, -0.002, -0.005],
            [-0.002,  0.148,  0.052, ...,  0.005, -0.006, -0.002],
            [   nan,    nan,    nan, ...,  0.002, -0.003, -0.005],
            ...,
            [   nan,    nan,    nan, ..., -0.014, -0.005, -0.021],
            [   nan,    nan,    nan, ..., -0.053, -0.025, -0.06 ],
            [   nan,    nan,    nan, ..., -0.011, -0.009, -0.02 ]],
    
           [[   nan,    nan,    nan, ..., -0.016, -0.005,  0.002],
            [   nan,    nan,    nan, ...,  0.011, -0.012,  0.006],
            [   nan,    nan,    nan, ...,  0.011, -0.002, -0.025],
            ...,
            [   nan,    nan,    nan, ..., -0.037, -0.043, -0.028],
            [   nan,    nan,    nan, ..., -0.044,  0.061, -0.001],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]]],
          dtype=float32)

  • amp
    (beam, range, time)
    uint8
    84 85 85 85 85 ... 148 154 161 147

    units :

    counts

    array([[[ 84,  85,  85, ..., 104, 105, 109],
            [ 49,  50,  50, ..., 101, 104, 105],
            [ 37,  37,  37, ...,  96,  98,  97],
            ...,
            [ 36,  36,  36, ..., 181, 181, 178],
            [ 36,  36,  36, ..., 161, 167, 163],
            [ 36,  36,  36, ..., 148, 158, 152]],
    
           [[111, 111, 111, ..., 104, 104, 106],
            [ 75,  75,  75, ..., 103, 105, 104],
            [ 44,  44,  44, ...,  96,  97,  98],
            ...,
            [ 34,  34,  34, ..., 209, 208, 209],
            [ 34,  34,  34, ..., 195, 195, 193],
            [ 34,  34,  34, ..., 168, 173, 164]],
    
           [[112, 112, 112, ..., 104, 109, 111],
            [ 75,  75,  75, ..., 103, 107, 107],
            [ 49,  49,  49, ...,  98, 100, 102],
            ...,
            [ 40,  40,  40, ..., 213, 214, 213],
            [ 40,  40,  40, ..., 188, 191, 184],
            [ 40,  40,  40, ..., 168, 172, 160]],
    
           [[110, 110, 110, ..., 105, 105, 108],
            [ 70,  70,  70, ..., 101, 107, 105],
            [ 40,  40,  40, ...,  97, 100,  99],
            ...,
            [ 35,  35,  35, ..., 198, 194, 193],
            [ 35,  35,  35, ..., 170, 171, 164],
            [ 35,  35,  35, ..., 154, 161, 147]]], dtype=uint8)

  • corr
    (beam, range, time)
    uint8
    68 67 67 66 65 ... 82 77 82 90 82

    units :

    counts

    array([[[ 68,  67,  67, ..., 119, 118, 117],
            [ 41,  42,  41, ..., 120, 123, 121],
            [ 23,  23,  21, ..., 126, 121, 124],
            ...,
            [ 20,  20,  20, ..., 133, 129, 131],
            [ 18,  19,  19, ...,  62,  69,  59],
            [  0,   0,   0, ...,  77,  84,  82]],
    
           [[ 62,  63,  62, ..., 123, 119, 116],
            [ 81,  77,  75, ..., 116, 121, 121],
            [ 44,  45,  45, ..., 125, 129, 132],
            ...,
            [ 19,  19,  19, ..., 104, 104, 100],
            [ 19,  19,  20, ...,  69,  79,  68],
            [ 19,  18,  18, ...,   0,   0,   0]],
    
           [[ 35,  35,  34, ..., 116, 119, 119],
            [ 37,  44,  47, ..., 120, 123, 119],
            [ 63,  62,  62, ..., 125, 123, 125],
            ...,
            [ 20,  19,  20, ..., 136, 135, 135],
            [ 20,  20,  20, ...,  77,  80,  69],
            [ 20,  19,  19, ...,  84,  88,  80]],
    
           [[ 63,  58,  60, ..., 122, 115, 118],
            [ 68,  70,  70, ..., 118, 122, 120],
            [ 35,  39,  39, ..., 119, 129, 123],
            ...,
            [ 19,  18,  18, ...,  91,  92,  84],
            [ 19,  17,  18, ...,  74,  77,  70],
            [ 18,  18,  18, ...,  82,  90,  82]]], dtype=uint8)

  • prcnt_gd
    (beam, range, time)
    uint8
    28 12 11 21 14 2 0 ... 0 0 0 0 0 0

    units :

    %

    array([[[ 28,  12,  11, ...,   3,   4,   0],
            [  0,   0,   1, ...,   3,   0,   0],
            [  0,   0,   0, ...,   3,   1,   0],
            ...,
            [  0,   0,   0, ...,  18,  15,  26],
            [  0,   0,   0, ...,  37,  38,  36],
            [  0,   0,   0, ...,  49,  67,  57]],
    
           [[  0,   0,   0, ...,   0,   0,   0],
            [  0,   0,   0, ...,   0,   0,   0],
            [  0,   0,   0, ...,   0,   0,   0],
            ...,
            [  0,   0,   0, ...,   1,   0,   1],
            [  0,   0,   0, ...,   0,   0,   0],
            [  0,   0,   0, ...,   0,   0,   0]],
    
           [[ 71,  87,  88, ...,   0,   0,   0],
            [ 99,  99,  98, ...,   0,   1,   0],
            [100, 100, 100, ...,   2,   0,   0],
            ...,
            [100, 100, 100, ...,   1,   0,   1],
            [100, 100, 100, ...,  44,  26,  54],
            [100, 100, 100, ...,  50,  32,  42]],
    
           [[  0,   0,   0, ...,  95,  95,  99],
            [  0,   0,   0, ...,  96,  97,  98],
            [  0,   0,   0, ...,  93,  98,  99],
            ...,
            [  0,   0,   0, ...,  78,  82,  70],
            [  0,   0,   0, ...,  18,  34,   9],
            [  0,   0,   0, ...,   0,   0,   0]]], dtype=uint8)

  • beam2inst_orientmat
    (x, x*)
    float64
    1.462 -1.462 0.0 ... -1.034 -1.034

    array([[ 1.4619022 , -1.4619022 ,  0.        ,  0.        ],
           [ 0.        ,  0.        , -1.4619022 ,  1.4619022 ],
           [ 0.26604444,  0.26604444,  0.26604444,  0.26604444],
           [ 1.03372096,  1.03372096, -1.03372096, -1.03372096]])

  • orientmat
    (earth, inst, time)
    float64
    0.3861 0.3858 ... -0.9983 -0.9983

    array([[[ 0.38607788,  0.38575819,  0.3857533 , ...,  0.79047287,
              0.79047287,  0.79047287],
            [ 0.92143434,  0.9215765 ,  0.92158675, ...,  0.610044  ,
              0.610044  ,  0.610044  ],
            [ 0.04361883,  0.04344447,  0.04327016, ..., -0.05476265,
             -0.05476265, -0.05476265]],
    
           [[ 0.92238986,  0.92252469,  0.92252553, ...,  0.60985976,
              0.60985976,  0.60985976],
            [-0.38622236, -0.38590029, -0.38590062, ..., -0.79220515,
             -0.79220515, -0.79220515],
            [-0.00540534, -0.00540539, -0.00523106, ..., -0.02195638,
             -0.02195638, -0.02195638]],
    
           [[ 0.0118659 ,  0.01178376,  0.01187711, ..., -0.05677761,
             -0.05677761, -0.05677761],
            [ 0.04232045,  0.04216377,  0.04193572, ..., -0.01604161,
             -0.01604161, -0.01604161],
            [-0.99903363, -0.9990412 , -0.99904966, ..., -0.99825799,
             -0.99825799, -0.99825799]]])

• Attributes: (39)

  inst_make :

  TRDI

  inst_model :

  Workhorse

  inst_type :

  ADCP

  rotate_vars :

  ['vel']

  has_imu :

  0

  name :

  wh-adcp

  sourceprog :

  instrument

  prog_ver :

  16.15

  beam_angle :

  20

  freq :

  300

  beam_pattern :

  convex

  orientation :

  up

  n_beams :

  4

  n_cells :

  40

  pings_per_ensemble :

  220

  cell_size :

  1.0

  blank :

  1.0

  prof_mode :

  1

  corr_threshold :

  64

  prof_codereps :

  3

  min_pgood :

  0

  evel_threshold :

  2000

  sec_between_ping_groups :

  1.0

  coord_sys :

  earth

  use_pitchroll :

  yes

  use_3beam :

  yes

  bin_mapping :

  yes

  xducer_misalign_deg :

  0.0

  magnetic_var_deg :

  0.0

  sensors_src :

  01111101

  sensors_avail :

  00011101

  bin1_dist_m :

  2.49

  xmit_pulse :

  1.47

  water_ref_cells :

  [1, 5]

  fls_target_threshold :

  50

  xmit_lag_m :

  0.5

  fs :

  0.004545454545454545

  vel_gps_corrected :

  0

  h_deploy :

  0.6

%matplotlib notebook
fig, ax = plt.subplots(figsize=(12,4))
ds['heading'].plot()
ax.grid(True);

Zoom in and find where heading is stable (in our case it is after 2021-8-16 14:30) , the same holds for end

%matplotlib inline
fig, ax = plt.subplots(figsize=(12,4))
ds['heading'].sel(time=slice('2021-8-16 14:30:00','2021-12-09 14:00:00')).plot()
ax.grid(True);

Cut for the time range

myds = ds.sel(range=slice(0,39), time=slice('2021-8-16 14:30:00','2021-12-09 14:00:00')).squeeze()

api.clean.set_range_offset(myds, 0.6) # instrument was 0.6m above the bottom (this is based on the mooring setup)
adcp_depth = 39.0  # based on notes from the deployment, we need this if there is no pressure sensor on ADCP
adcp_lon = 16.3882 
adcp_lat = 43.5190

dlfn.set_declination(myds, declin=4.5, inplace=True) # 4 deg 30 sec East for Split

QC procedure

1. start with percentage of good solutions for 3 and 4 bins (>50%)
2. continue with correlation (bin3 or bin4 > 110)
3. amplitude of the signal (detect surface)

Percentage of good 3 & 4 bin solution mask

pg = myds['prcnt_gd'].values # % of 4 bin solutions

np.shape(pg)

(4, 36, 11039)

pg1 = pg[0,:,:] # % of 3 bin good solutions
pg4 = pg[3,:,:] # % of 3 bin good solutions
pg_test = pg1 + pg4
drange = myds['range'].values
t = dlfn.time.dt642date(myds.time)

%matplotlib inline
fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, pg_test, cmap='gist_ncar', vmin=0, vmax=100)
ax.set_title('Percentage of good solutions for bin 3 & 4')
fig.colorbar(c, ax=ax)
ax.grid(True)

# if pg1 + Pg4 < 50% that is definitely bad bin/data
mask_pg = np.ma.masked_invalid(pg_test>50)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask_pg)
ax.set_title('Percent good mask')
fig.colorbar(c, ax=ax)
ax.grid(True)

Bin correlation mask

# correlation btw beams, we can get surface reflection
myds['corr'].sel(beam=2).squeeze().plot(vmin = 0, vmax =110);

corr_3 = myds['corr'].sel(beam=2).squeeze().values # 3 bin 
corr_4 = myds['corr'].sel(beam=3).squeeze().values # 4 bin

# correlation QC if corr_3 or corr_4 < 110 
mask_corr = np.ma.masked_invalid(corr_3>110)
mask_corr = np.ma.masked_invalid(corr_4>110)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask_corr)
ax.set_title('Correlation mask')
fig.colorbar(c, ax=ax)
ax.grid(True)

Echo intensity (amplitude) mask

# reflection amplitude, we can see surface reflection of high turbidity
myds['amp'].sel(beam=3).squeeze().plot(cmap='gist_ncar');

amp_1 = myds['amp'].sel(beam=1).squeeze().values # 1 bin 
amp_2 = myds['amp'].sel(beam=2).squeeze().values # 2 bin 
amp_3 = myds['amp'].sel(beam=3).squeeze().values # 3 bin
mask_amp = np.ma.masked_invalid(amp_1<180)
mask_amp = np.ma.masked_invalid(amp_2<180)
mask_amp = np.ma.masked_invalid(amp_3<180)
amp_mean=1./3*(amp_1+amp_2+amp_3)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange[:-1]-adcp_depth, np.diff(amp_mean, axis=0), vmin = 0, vmax=30)
ax.set_title('Amplitude (echo intensity) diff')
fig.colorbar(c, ax=ax)
ax.grid(True)

mask_amp[0:-1,:] = np.ma.masked_invalid(np.diff(amp_mean, axis=0) <30)  # if change > 30 counts btw successive bins

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask_amp, vmin=0, vmax = 1)
ax.set_title('Amplitude (echo intensity) mask')
fig.colorbar(c, ax=ax)
ax.grid(True)

mask = mask_pg * mask_corr * mask_amp # stack all together

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask, vmin=0, vmax = 1)
ax.set_title('Combined mask (PG + CORR + AMP)')
fig.colorbar(c, ax=ax)
ax.grid(True)

Velocity limits mask

vel_mag = myds.velds.U_mag
val_dir = myds.velds.U_dir
vel_east = myds.velds.u
vel_north = myds.velds.v
vel_up = myds.velds.w

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, vel_up, vmin=0, vmax = 0.03)
ax.set_title('Upward velocity')
fig.colorbar(c, ax=ax)
ax.grid(True)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, vel_mag, vmin=0, vmax = 0.3)
ax.set_title('Horizontal velocity magnitude')
fig.colorbar(c, ax=ax)
ax.grid(True)

mask_vel = np.ma.masked_invalid(vel_up < 0.03)
mask_vel = np.ma.masked_invalid(vel_mag < 0.8)
mask_vel[0:-1,:] = np.ma.masked_invalid(np.diff(vel_mag, axis=0) <0.1) # if too large shear its bad
mask_vel[-2:,:] = 0 # last 2 bins are bad as too close to surface (40 bins in 39m water)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask_vel, vmin=0, vmax = 1)
ax.set_title('Velocity mask')
fig.colorbar(c, ax=ax)
ax.grid(True)

mask = mask * mask_vel # stack masks

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, mask, vmin=0, vmax = 1)
ax.set_title('Combined final mask')
fig.colorbar(c, ax=ax)
ax.grid(True)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, np.ma.masked_where(~mask, vel_east),cmap='bwr', vmin=-0.3, vmax = 0.3)
ax.set_title('East component')
fig.colorbar(c, ax=ax)
ax.grid(True)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, np.ma.masked_where(~mask, vel_north),cmap='bwr', vmin=-0.3, vmax = 0.3)
ax.set_title('North component')
fig.colorbar(c, ax=ax)
ax.grid(True)

# fill small gaps (limit for maxgap of 3 nans)
north = np.ma.masked_where(~mask, vel_north)
north = north.filled(np.nan)
dlfn.tools.misc.fillgaps(north, dim=0, maxgap = 3)
dlfn.tools.misc.fillgaps(north, dim=1, maxgap = 3)

east = np.ma.masked_where(~mask, vel_east)
east = east.filled(np.nan)
dlfn.tools.misc.fillgaps(east, dim=0, maxgap = 3)
dlfn.tools.misc.fillgaps(east, dim=1, maxgap = 3)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, north,cmap='bwr', vmin=-0.3, vmax = 0.3)
ax.set_title('North component')
fig.colorbar(c, ax=ax)
ax.grid(True)

fig, ax = plt.subplots(figsize=(13,5))
c = ax.pcolormesh(t, drange-adcp_depth, east, cmap='bwr', vmin=-0.3, vmax = 0.3)
ax.set_title('East component')
ax.set_xlim()
fig.colorbar(c, ax=ax)
ax.grid(True)

temperature =  myds['temp'].values

fig, ax = plt.subplots(figsize=(13,5))
ax.plot(t, temperature, color = 'k')
ax.set_title('Temperature at the bottom (38.4m)')
ax.set_xlim([t[0], t[-1]])
ax.grid(True)

# save that to netcdf file, no need for xarray as I have only 4 variables
outfile = 'ADCP_CAAT_01.nc'
!rm -f ADCP_CAAT_01.nc
from netCDF4 import Dataset, num2date, date2num
nf_out = Dataset(outfile, "w", clobber = True , format="NETCDF4_CLASSIC")
nf_out.title = "Analyzed data for ADCP 01 (lon = 16.3882, lat = 43.519, depth = 39m)"
nf_out.description = "ADCP deployed by IOR *Hrvoje Mihanovic* on 2021-08-16"
nf_out.author = "Analysis done by Ivica Janekovic"
nf_out.type = "netCDF4"
# dimensions
nf_out.createDimension('time', None)
nf_out.createDimension('bin', north.shape[0])
nf_out.createDimension('one',1)
# ADCP station depth
depth = nf_out.createVariable('depth', 'f', ('one'),  zlib=False)
depth.units = 'm'
depth.long_name = 'ADCP station depth'
depth.short_name = 'depth'
depth[:] = adcp_depth
# ADCP station longitude
lon = nf_out.createVariable('lon', 'd', ('one'),  zlib=False)
lon.units = 'degrees east'
lon.long_name = 'Longitude for ADCP location'
lon.short_name ='longitude'
lon[:] = adcp_lon
# ADCP station latitude
lat = nf_out.createVariable('lat', 'd', ('one'),  zlib=False)
lat.units = 'degrees north'
lat.long_name = 'Latitude for ADCP location'
lat.short_name ='latitude'
lat[:] = adcp_lat
# ADCP station time
time = nf_out.createVariable('time', 'd', ('time'),  zlib=False)
time.units = "days since 2000-01-01 00:00:00"
time.calendar = 'proleptic_gregorian'
time.long_name = 'ADCP time (UTC)'
time.standard_name = 'time'
time.coordinate_type = 'time'
time[:] = date2num(t, time.units)
# Temperature at the ADCP (bottom)
temp = nf_out.createVariable('temp', 'f', ('time'),  zlib=True, least_significant_digit=3, complevel=7)
temp.units = 'C'
temp.coordinates = 'time'
temp.long_name = 'Ocean temperature at the bottom (ADCP location)'
temp.short_name ='temp'
temp[:] = temperature
# Beam depths (from the free surface)
beam_depth = nf_out.createVariable('beam_depth', 'f', ('bin'),  zlib=False)
beam_depth.units = "m"
beam_depth.long_name = 'ADCP beam depths (from the free surace), negative values in [m]'
beam_depth.standard_name = 'depth'
beam_depth.coordinate_type = 'bin'
beam_depth[:] =  drange - adcp_depth
# Eastward current
eastward_current = nf_out.createVariable('eastward_current', 'f', ('bin','time'), fill_value=1.e+20,  
                                         zlib=True, least_significant_digit=3, complevel=7)
eastward_current.units = "ms^-1"
eastward_current.long_name = 'ADCP eastward component of ocean current'
eastward_current.standard_name = 'east_current'
eastward_current.coordinate_type = 'time, bin'
eastward_current[:] = east
# Northward current
northward_current = nf_out.createVariable('northward_current', 'f', ('bin','time'), fill_value = 1.e+20,  
                                         zlib=True, least_significant_digit=3, complevel=7)
northward_current.units = "ms^-1"
northward_current.long_name = 'ADCP northward component of ocean current'
northward_current.standard_name = 'north_current'
northward_current.coordinate_type = 'time, bin'
northward_current[:] = north
nf_out.close()