diff --git a/impdar/bin/impplot.py b/impdar/bin/impplot.py
index 0ca46e94..c433a3ec 100644
--- a/impdar/bin/impplot.py
+++ b/impdar/bin/impplot.py
@@ -24,10 +24,19 @@ def _get_args():
                                        plot_radargram,
                                        defname='radargram',
                                        xd=True,
-                                       yd=True)
-    rg_parser.add_argument('-picks', action='store_true', help='Plot picks')
-    rg_parser.add_argument('-clims', nargs=2, type=float, help='Color limits')
-    rg_parser.add_argument('-flatten_layer', type=int, default=None, help='Distort plot so this layer is flat')
+                                       yd=True,
+                                       dualy=True)
+    rg_parser.add_argument('-picks',
+                           action='store_true',
+                           help='Plot picks')
+    rg_parser.add_argument('-clims',
+                           nargs=2,
+                           type=float,
+                           help='Color limits')
+    rg_parser.add_argument('-flatten_layer',
+                           type=int,
+                           default=None,
+                           help='Distort plot so this layer is flat')
     rg_parser.add_argument('-cmap',
                            type=str,
                            default='gray',
@@ -37,17 +46,13 @@ def _get_args():
                            'ft',
                            'Plot ft',
                            plot_ft,
-                           defname='spec',
-                           xd=True,
-                           yd=True)
+                           defname='spec')
 
     _add_simple_procparser(subparsers,
                            'hft',
                            'Plot ft',
                            plot_hft,
-                           defname='spec',
-                           xd=True,
-                           yd=True)
+                           defname='spec')
 
     trace_parser = _add_simple_procparser(subparsers,
                                           'traces',
@@ -55,19 +60,20 @@ def _get_args():
                                           plot_traces,
                                           defname='traces',
                                           xd=False,
-                                          yd=True)
+                                          yd=True,
+                                          dualy=True)
     trace_parser.add_argument('t_start',
                               type=int,
                               help='Starting trace number')
-    trace_parser.add_argument('t_end', type=int, help='Ending trace number')
+    trace_parser.add_argument('t_end',
+                              type=int,
+                              help='Ending trace number')
 
     power_parser = _add_simple_procparser(subparsers,
                                           'power',
                                           'Plot power on a layer',
                                           plot_power,
                                           defname='power',
-                                          xd=False,
-                                          yd=False,
                                           other_ftypes=False)
     power_parser.add_argument('layer',
                               type=int,
@@ -78,8 +84,6 @@ def _get_args():
                                          'Plot spectrogram for all traces',
                                          plot_spectrogram,
                                          defname='spectrogram',
-                                         xd=False,
-                                         yd=False,
                                          other_ftypes=False)
     spec_parser.add_argument('freq_lower',
                              type=float,
@@ -91,10 +95,10 @@ def _get_args():
 
 
 def _add_simple_procparser(subparsers, name, helpstr, func, defname='proc',
-                           xd=False, yd=False, other_ftypes=True):
+                           xd=False, yd=False, dualy=False, other_ftypes=True):
     """Add a simple subparser."""
     parser = _add_procparser(subparsers, name, helpstr, func, defname=defname)
-    _add_def_args(parser, xd=xd, yd=yd)
+    _add_def_args(parser, xd=xd, yd=yd, dualy=dualy)
     return parser
 
 
@@ -105,7 +109,7 @@ def _add_procparser(subparsers, name, helpstr, func, defname='proc'):
     return parser
 
 
-def _add_def_args(parser, xd=False, yd=False, other_ftypes=True):
+def _add_def_args(parser, xd=False, yd=False, dualy=False, other_ftypes=True):
     """Set up common arguments for different types of commands."""
     parser.add_argument('fns',
                         type=str,
@@ -136,6 +140,11 @@ def _add_def_args(parser, xd=False, yd=False, other_ftypes=True):
                             action='store_true',
                             help='Plot the depth rather than travel time')
 
+    if dualy:
+        parser.add_argument('-dualy',
+                            action='store_true',
+                            help='Primary y axis is TWTT, secondary is depth')
+
     if other_ftypes:
         parser.add_argument('--in_fmt', type=str,
                             help='Type of file',
@@ -145,11 +154,11 @@ def _add_def_args(parser, xd=False, yd=False, other_ftypes=True):
 
 def plot_radargram(fns=None, s=False, o=None, xd=False, yd=False, o_fmt='png',
                    dpi=300, in_fmt='mat', picks=False, clims=None, cmap='gray',
-                   flatten_layer=None, **kwargs):
+                   flatten_layer=None, dualy=False, **kwargs):
     """Plot data as a radio echogram."""
     plot.plot(fns, xd=xd, yd=yd, s=s, o=o, ftype=o_fmt, dpi=dpi,
               filetype=in_fmt, pick_colors=picks, cmap=cmap, clims=clims,
-              flatten_layer=flatten_layer)
+              flatten_layer=flatten_layer, dualy=dualy)
 
 
 def plot_ft(fns=None, s=False, o=None, xd=False, yd=False, o_fmt='png',
@@ -181,11 +190,11 @@ def plot_power(fns=None, layer=None, s=False, o=None, o_fmt='png',
               filetype=in_fmt)
 
 
-def plot_traces(fns=None, t_start=None, t_end=None, yd=False, s=False, o=None,
-                o_fmt='png', dpi=300, in_fmt='mat', **kwargs):
+def plot_traces(fns=None, t_start=None, t_end=None, yd=False, dualy=False,
+                s=False, o=None, o_fmt='png', dpi=300, in_fmt='mat', **kwargs):
     """Plot traces in terms of amplitude vs some vertical variable."""
     plot.plot(fns, tr=(t_start, t_end), yd=yd, s=s, o=o, ftype=o_fmt, dpi=dpi,
-              filetype=in_fmt)
+              dualy=dualy, filetype=in_fmt)
 
 
 def plot_spectrogram(fns=None, freq_lower=None, freq_upper=None, window=None,
diff --git a/impdar/lib/RadarData/_RadarDataProcessing.py b/impdar/lib/RadarData/_RadarDataProcessing.py
index aeae0791..3ad19528 100644
--- a/impdar/lib/RadarData/_RadarDataProcessing.py
+++ b/impdar/lib/RadarData/_RadarDataProcessing.py
@@ -190,7 +190,7 @@ def optimize_moveout_depth(d_in, t, ant_sep, profile_depth, profile_u):
     profile_u: array
         velocity
     """
-    args = np.argwhere(profile_depth<d_in)
+    args = np.argwhere(profile_depth<=d_in)
     if len(args) == 0:
         raise ValueError('Profile not shallow enough. Extend to cover top')
     vels = profile_u[args][:,0]
diff --git a/impdar/lib/load/__init__.py b/impdar/lib/load/__init__.py
index f5fb0de1..13722fc9 100644
--- a/impdar/lib/load/__init__.py
+++ b/impdar/lib/load/__init__.py
@@ -12,8 +12,9 @@
 
 import os.path
 import numpy as np
-from . import load_gssi, load_pulse_ekko, load_gprMax, load_olaf, load_mcords, load_segy, load_UoA_mat, load_ramac, load_bsi
-from . import load_delores, load_osu, load_stomat
+from . import load_mcords  # needs to be imported first and alone due to opaque h5py/netcdf4 error
+from . import load_gssi, load_pulse_ekko, load_gprMax, load_olaf, load_segy, load_UoA_mat
+from . import load_delores, load_osu, load_stomat, load_ramac, load_bsi
 from ..RadarData import RadarData
 
 # This should be updated as new functionality arrives
diff --git a/impdar/lib/plot.py b/impdar/lib/plot.py
index 81c61282..1edab5be 100644
--- a/impdar/lib/plot.py
+++ b/impdar/lib/plot.py
@@ -18,10 +18,10 @@
                   '#882255', '#44AA99', '#999933', '#AA4499']
 
 def plot(fns, tr=None, s=False, ftype='png', dpi=300, xd=False, yd=False,
-         x_range=(0, -1), power=None, spectra=None, freq_limit=None,
-         window=None, scaling='spectrum', filetype='mat', pick_colors=None,
-         ft=False, hft=False, clims=None, cmap=plt.cm.gray, flatten_layer=None,
-         *args, **kwargs):
+         dualy=False, x_range=(0, -1), power=None, spectra=None,
+         freq_limit=None, window=None, scaling='spectrum', filetype='mat',
+         pick_colors=None, ft=False, hft=False, clims=None, cmap=plt.cm.gray,
+         flatten_layer=None, *args, **kwargs):
     """Wrap a number of plot types.
 
     This should really only be used by the exectuables.
@@ -52,7 +52,11 @@ def plot(fns, tr=None, s=False, ftype='png', dpi=300, xd=False, yd=False,
     else:
         xdat = 'tnum'
     if yd:
+        if dualy:
+            raise ValueError('Only one of yd and dualy can be true')
         ydat = 'depth'
+    elif dualy:
+        ydat = 'dual'
     else:
         ydat = 'twtt'
 
@@ -164,7 +168,8 @@ def plot_radargram(dat, xdat='tnum', ydat='twtt', x_range=(0, -1),
 
     if y_range is None:
         y_range = (0, -1)
-    elif y_range[-1] == -1:
+
+    if y_range[-1] == -1:
         y_range = (y_range[0], dat.data.shape[0])
 
     if dat.data.dtype in [np.complex128]:
@@ -197,6 +202,8 @@ def norm(x):
     else:
         ax.invert_yaxis()
         if ydat == 'twtt':
+            # we have a chance that there are NaNs after NMO correction...
+            y_range = (max(y_range[0], np.min(np.where(~np.isnan(dat.travel_time))[0])), y_range[1])
             yd = dat.travel_time[y_range[0]:y_range[-1]]
             ax.set_ylabel('Two way travel time (usec)')
         elif ydat == 'depth':
@@ -206,9 +213,23 @@ def norm(x):
                 yd = dat.travel_time[y_range[0]:y_range[-1]] / 2.0 * (
                     1.69e8 * 1.0e-6)
             ax.set_ylabel('Depth (m)')
+        elif ydat == 'dual':
+            # we have a chance that there are NaNs after NMO correction...
+            y_range = (max(y_range[0], np.min(np.where(~np.isnan(dat.travel_time))[0])), y_range[1])
+
+            yd = dat.travel_time[y_range[0]:y_range[-1]]
+            ax.set_ylabel('Two way travel time (usec)')
+            ax2 = ax.twinx()
+            if dat.nmo_depth is not None:
+                yd2 = dat.nmo_depth[y_range[0]:y_range[-1]]
+            else:
+                yd2 = dat.travel_time[y_range[0]:y_range[-1]] / 2.0 * (
+                    1.69e8 * 1.0e-6)
+            ax2.set_ylabel('Approximate depth (m)')
+            ax2.set_ylim(yd2[-1], yd2[0])
         else:
             raise ValueError('Unrecognized ydat, choices are elev, twtt, \
-                             or depth')
+                             depth, or dual')
 
     if xdat == 'tnum':
         xd = np.arange(int(dat.tnum))[x_range[0]:x_range[-1]]
@@ -379,7 +400,7 @@ def plot_traces(dat, tr, ydat='twtt', fig=None, ax=None, linewidth=1.0,
     elif tr[0] == tr[1]:
         tr = (tr[0], tr[0] + 1)
 
-    if ydat not in ['twtt', 'depth']:
+    if ydat not in ['twtt', 'depth', 'dual']:
         raise ValueError('y axis choices are twtt or depth')
     if fig is not None:
         if ax is None:
@@ -401,6 +422,19 @@ def plot_traces(dat, tr, ydat='twtt', fig=None, ax=None, linewidth=1.0,
         else:
             yd = dat.nmo_depth
         ax.set_ylabel('Depth (m)')
+    elif ydat == 'dual':
+        # we have a chance that there are NaNs after NMO correction...
+        yd = dat.travel_time
+        ax.set_ylabel('Two way travel time (usec)')
+        ax2 = ax.twinx()
+        if dat.nmo_depth is not None:
+            yd2 = dat.nmo_depth
+        else:
+            yd2 = dat.travel_time / 2.0 * (1.69e8 * 1.0e-6)
+        ax2.set_ylabel('Approximate depth (m)')
+        ax2.set_ylim(yd2[-1], yd2[0])
+    else:
+        raise ValueError('Unrecognized y scale')
 
     for j in range(*tr):
         ax.plot(dat.data[:, j], yd, linewidth=linewidth, linestyle=linestyle)
diff --git a/impdar/tests/input_data/rho_profile.txt b/impdar/tests/input_data/rho_profile.txt
index 10cc2f70..2b53167f 100644
--- a/impdar/tests/input_data/rho_profile.txt
+++ b/impdar/tests/input_data/rho_profile.txt
@@ -1,5 +1,5 @@
-800,0
-900,50
-910,51
-910,100
-910,10000
+0,800
+50,900
+51,910
+100,910
+10000,910
diff --git a/impdar/tests/test_RadarData.py b/impdar/tests/test_RadarData.py
index c985c9f8..721989d8 100644
--- a/impdar/tests/test_RadarData.py
+++ b/impdar/tests/test_RadarData.py
@@ -185,10 +185,8 @@ def test_rangegain(self):
         self.data.rangegain(1.0)
         self.assertTrue(self.data.flags.rgain)
 
-    @patch('impdar.lib.RadarData._RadarDataProcessing.optimize_moveout_depth', returns=[1000.])
-    def test_NMO(self, mock_omd):
+    def test_NMO(self):
         # If velocity is 2
-
         self.data.nmo(0., uice=2.0, uair=2.0)
         self.assertTrue(np.allclose(self.data.travel_time * 1.0e-6, self.data.nmo_depth))
 
diff --git a/impdar/tests/test_plot.py b/impdar/tests/test_plot.py
index 3e4fce4d..b14e6dfa 100644
--- a/impdar/tests/test_plot.py
+++ b/impdar/tests/test_plot.py
@@ -138,10 +138,12 @@ def test_plot_traces(self, mock_show):
 
         # no nmo
         plot.plot_traces(dat, 0, ydat='depth', fig=fig, ax=ax)
+        plot.plot_traces(dat, 0, ydat='dual', fig=fig, ax=ax)
 
         # with nmo
-        dat.nmo_depth = np.arange(10)
+        dat.nmo_depth = np.linspace(0, 10, dat.travel_time.shape[0])
         plot.plot_traces(dat, 0, ydat='depth', fig=fig, ax=ax)
+        plot.plot_traces(dat, 0, ydat='dual', fig=fig, ax=ax)
         with self.assertRaises(ValueError):
             plot.plot_traces(dat, 0, ydat='dum', fig=fig, ax=ax)
 
@@ -219,8 +221,18 @@ def test_plot_radargram(self, mock_show):
         with self.assertRaises(ValueError):
             plot.plot_radargram(dat, xdat='dummy', fig=fig, ax=ax)
 
+        # Varying ydata
+        dat.nmo_depth = None
         plot.plot_radargram(dat, y_range=None, fig=fig, ax=ax)
         plot.plot_radargram(dat, ydat='depth', fig=fig, ax=ax)
+        plot.plot_radargram(dat, ydat='dual', fig=fig, ax=ax)
+
+        # with nmo defined, these two are different
+        dat.nmo_depth = np.linspace(0, 100, dat.travel_time.shape[0])
+        plot.plot_radargram(dat, ydat='depth', fig=fig, ax=ax)
+        plot.plot_radargram(dat, ydat='dual', fig=fig, ax=ax)
+
+        dat = NoInitRadarData(big=True)
         with self.assertRaises(ValueError):
             plot.plot_radargram(dat, ydat='dummy', fig=fig, ax=ax)