Skip to content

Commit

Permalink
Dual axes (#14)
Browse files Browse the repository at this point in the history 
* Add support for dual axes (TWTT and depth)

* Fix very opaque error: import order of h5py and netCDF4 matters.
See comment on Unidata/netcdf4-python#653
by arthur-e

* Fix syntax error

* Fix bug in variable-rho NMO exposed by update to scipy
  • Loading branch information
@dlilien
dlilien committed Jul 2, 2020
1 parent 9467068 commit 00c6e86
Show file tree
Hide file tree
Showing 7 changed files with 97 additions and 43 deletions.
57 changes: 33 additions & 24 deletions impdar/bin/impplot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand All @@ -37,37 +46,34 @@ 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',
'Plot traces vs depth',
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,
Expand All @@ -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,
Expand All @@ -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


Expand All @@ -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,
Expand Down Expand Up @@ -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',
Expand All @@ -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',
Expand Down Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion impdar/lib/RadarData/_RadarDataProcessing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
5 changes: 3 additions & 2 deletions impdar/lib/load/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
48 changes: 41 additions & 7 deletions impdar/lib/plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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'

Expand Down Expand Up @@ -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]:
Expand Down Expand Up @@ -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':
Expand All @@ -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]]
Expand Down Expand Up @@ -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:
Expand All @@ -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)
Expand Down
10 changes: 5 additions & 5 deletions impdar/tests/input_data/rho_profile.txt
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
800,0
900,50
910,51
910,100
910,10000
0,800
50,900
51,910
100,910
10000,910
4 changes: 1 addition & 3 deletions impdar/tests/test_RadarData.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Expand Down
14 changes: 13 additions & 1 deletion impdar/tests/test_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down

0 comments on commit 00c6e86

Please sign in to comment.