你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 04:56:51 +08:00
Updated so they can be imported or run as scripts from the command line.
这个提交包含在:
@@ -25,188 +25,201 @@ import matplotlib.gridspec as gridspec
|
|||||||
from gprMax.exceptions import CmdInputError
|
from gprMax.exceptions import CmdInputError
|
||||||
from gprMax.receivers import Rx
|
from gprMax.receivers import Rx
|
||||||
|
|
||||||
"""Plots electric and magnetic fields and currents from all receiver points in the given output file. Each receiver point is plotted in a new figure window."""
|
|
||||||
|
|
||||||
# Parse command line arguments
|
def make_plot(filename, outputs=Rx.availableoutputs, fft=False):
|
||||||
parser = argparse.ArgumentParser(description='Plots electric and magnetic fields and currents from all receiver points in the given output file. Each receiver point is plotted in a new figure window.', usage='cd gprMax; python -m tools.plot_Ascan outputfile')
|
"""Plots electric and magnetic fields and currents from all receiver points in the given output file. Each receiver point is plotted in a new figure window.
|
||||||
parser.add_argument('outputfile', help='name of output file including path')
|
|
||||||
parser.add_argument('--outputs', help='outputs to be plotted', choices='Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, Iz', default=Rx.availableoutputs, nargs='+')
|
Args:
|
||||||
parser.add_argument('-fft', action='store_true', default=False, help='plot FFT (single output must be specified)')
|
filename (string): Filename (including path) of output file.
|
||||||
args = parser.parse_args()
|
outputs (list): List of field/current components to plot.
|
||||||
|
fft (boolean): Plot FFT switch.
|
||||||
# Open output file and read some attributes
|
"""
|
||||||
f = h5py.File(args.outputfile, 'r')
|
|
||||||
nrx = f.attrs['nrx']
|
|
||||||
dt = f.attrs['dt']
|
|
||||||
iterations = f.attrs['Iterations']
|
|
||||||
time = np.linspace(0, 1, iterations)
|
|
||||||
time *= (iterations * dt)
|
|
||||||
|
|
||||||
# Check there are any receivers
|
|
||||||
if nrx == 0:
|
|
||||||
raise CmdInputError('No receivers found in {}'.format(args.outputfile))
|
|
||||||
|
|
||||||
# Check for single output component when doing a FFT
|
|
||||||
if args.fft:
|
|
||||||
if not len(args.outputs) == 1:
|
|
||||||
raise CmdInputError('A single output must be specified when using the -fft option')
|
|
||||||
|
|
||||||
# New plot for each receiver
|
|
||||||
for rx in range(1, nrx + 1):
|
|
||||||
path = '/rxs/rx' + str(rx) + '/'
|
|
||||||
availableoutputs = list(f[path].keys())
|
|
||||||
|
|
||||||
# If only a single output is required, create one subplot
|
# Open output file and read some attributes
|
||||||
if len(args.outputs) == 1:
|
f = h5py.File(filename, 'r')
|
||||||
|
nrx = f.attrs['nrx']
|
||||||
|
dt = f.attrs['dt']
|
||||||
|
iterations = f.attrs['Iterations']
|
||||||
|
time = np.linspace(0, 1, iterations)
|
||||||
|
time *= (iterations * dt)
|
||||||
|
|
||||||
|
# Check there are any receivers
|
||||||
|
if nrx == 0:
|
||||||
|
raise CmdInputError('No receivers found in {}'.format(filename))
|
||||||
|
|
||||||
|
# Check for single output component when doing a FFT
|
||||||
|
if fft:
|
||||||
|
if not len(outputs) == 1:
|
||||||
|
raise CmdInputError('A single output must be specified when using the -fft option')
|
||||||
|
|
||||||
|
# New plot for each receiver
|
||||||
|
for rx in range(1, nrx + 1):
|
||||||
|
path = '/rxs/rx' + str(rx) + '/'
|
||||||
|
availableoutputs = list(f[path].keys())
|
||||||
|
|
||||||
# Check for polarity of output and if requested output is in file
|
# If only a single output is required, create one subplot
|
||||||
if args.outputs[0][0] == 'm':
|
if len(outputs) == 1:
|
||||||
polarity = -1
|
|
||||||
outputtext = '-' + args.outputs[0][1:]
|
|
||||||
output = args.outputs[0][1:]
|
|
||||||
else:
|
|
||||||
polarity = 1
|
|
||||||
outputtext = args.outputs[0]
|
|
||||||
output = args.outputs[0]
|
|
||||||
|
|
||||||
if output not in availableoutputs:
|
|
||||||
raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(output, ', '.join(availableoutputs)))
|
|
||||||
|
|
||||||
outputdata = f[path + output][:] * polarity
|
|
||||||
|
|
||||||
# Plotting if FFT required
|
|
||||||
if args.fft:
|
|
||||||
# Calculate magnitude of frequency spectra of waveform
|
|
||||||
power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)
|
|
||||||
freqs = np.fft.fftfreq(power.size, d=dt)
|
|
||||||
|
|
||||||
# Shift powers so that frequency with maximum power is at zero decibels
|
|
||||||
power -= np.amax(power)
|
|
||||||
|
|
||||||
# Set plotting range to -60dB from maximum power
|
|
||||||
pltrange = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
|
|
||||||
# To a maximum frequency
|
|
||||||
#pltrange = np.where(freqs > 2e9)[0][0]
|
|
||||||
pltrange = np.s_[0:pltrange]
|
|
||||||
|
|
||||||
# Plot time history of output component
|
|
||||||
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
|
||||||
line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext)
|
|
||||||
ax1.set_xlabel('Time [s]')
|
|
||||||
ax1.set_ylabel(outputtext + ' field strength [V/m]')
|
|
||||||
ax1.set_xlim([0, np.amax(time)])
|
|
||||||
ax1.grid()
|
|
||||||
|
|
||||||
# Plot frequency spectra
|
|
||||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.')
|
|
||||||
plt.setp(baseline, 'linewidth', 0)
|
|
||||||
plt.setp(stemlines, 'color', 'r')
|
|
||||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
|
||||||
line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
|
||||||
ax2.set_xlabel('Frequency [Hz]')
|
|
||||||
ax2.set_ylabel('Power [dB]')
|
|
||||||
ax2.grid()
|
|
||||||
|
|
||||||
# Change colours and labels for magnetic field components or currents
|
|
||||||
if 'H' in args.outputs[0]:
|
|
||||||
plt.setp(line1, color='g')
|
|
||||||
plt.setp(line2, color='g')
|
|
||||||
plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]')
|
|
||||||
plt.setp(stemlines, 'color', 'g')
|
|
||||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
|
||||||
elif 'I' in args.outputs[0]:
|
|
||||||
plt.setp(line1, color='b')
|
|
||||||
plt.setp(line2, color='b')
|
|
||||||
plt.setp(ax1, ylabel=outputtext + ' current [A]')
|
|
||||||
plt.setp(stemlines, 'color', 'b')
|
|
||||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
|
||||||
|
|
||||||
plt.show()
|
|
||||||
|
|
||||||
# Plotting if no FFT required
|
|
||||||
else:
|
|
||||||
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]', ylabel=outputtext + ' field strength [V/m]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
|
||||||
line = ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
|
||||||
ax.set_xlim([0, np.amax(time)])
|
|
||||||
#ax.set_ylim([-15, 20])
|
|
||||||
ax.grid()
|
|
||||||
|
|
||||||
if 'H' in output:
|
|
||||||
plt.setp(line, color='g')
|
|
||||||
plt.setp(ax, ylabel=outputtext + ', field strength [A/m]')
|
|
||||||
elif 'I' in output:
|
|
||||||
plt.setp(line, color='b')
|
|
||||||
plt.setp(ax, ylabel=outputtext + ', current [A]')
|
|
||||||
|
|
||||||
# If multiple outputs required, create all nine subplots and populate only the specified ones
|
|
||||||
else:
|
|
||||||
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
|
||||||
gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3)
|
|
||||||
for output in args.outputs:
|
|
||||||
|
|
||||||
# Check for polarity of output and if requested output is in file
|
# Check for polarity of output and if requested output is in file
|
||||||
if output[0] == 'm':
|
if outputs[0][0] == 'm':
|
||||||
polarity = -1
|
polarity = -1
|
||||||
outputtext = '-' + output[1:]
|
outputtext = '-' + outputs[0][1:]
|
||||||
output = output[1:]
|
output = outputs[0][1:]
|
||||||
else:
|
else:
|
||||||
polarity = 1
|
polarity = 1
|
||||||
outputtext = output
|
outputtext = outputs[0]
|
||||||
|
output = outputs[0]
|
||||||
|
|
||||||
# Check if requested output is in file
|
|
||||||
if output not in availableoutputs:
|
if output not in availableoutputs:
|
||||||
raise CmdInputError('Output(s) requested to plot: {}, but available output(s) for receiver {} in the file: {}'.format(', '.join(args.outputs), rx, ', '.join(availableoutputs)))
|
raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(output, ', '.join(availableoutputs)))
|
||||||
|
|
||||||
outputdata = f[path + output][:] * polarity
|
outputdata = f[path + output][:] * polarity
|
||||||
|
|
||||||
if output == 'Ex':
|
|
||||||
ax = plt.subplot(gs[0, 0])
|
|
||||||
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
|
||||||
#ax.set_ylim([-15, 20])
|
|
||||||
elif output == 'Ey':
|
|
||||||
ax = plt.subplot(gs[1, 0])
|
|
||||||
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
|
||||||
#ax.set_ylim([-15, 20])
|
|
||||||
elif output == 'Ez':
|
|
||||||
ax = plt.subplot(gs[2, 0])
|
|
||||||
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
|
||||||
#ax.set_ylim([-15, 20])
|
|
||||||
elif output == 'Hx':
|
|
||||||
ax = plt.subplot(gs[0, 1])
|
|
||||||
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
|
||||||
#ax.set_ylim([-0.03, 0.03])
|
|
||||||
elif output == 'Hy':
|
|
||||||
ax = plt.subplot(gs[1, 1])
|
|
||||||
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
|
||||||
#ax.set_ylim([-0.03, 0.03])
|
|
||||||
elif output == 'Hz':
|
|
||||||
ax = plt.subplot(gs[2, 1])
|
|
||||||
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
|
||||||
#ax.set_ylim([-0.03, 0.03])
|
|
||||||
elif output == 'Ix':
|
|
||||||
ax = plt.subplot(gs[0, 2])
|
|
||||||
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', current [A]')
|
|
||||||
elif output == 'Iy':
|
|
||||||
ax = plt.subplot(gs[1, 2])
|
|
||||||
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', current [A]')
|
|
||||||
elif output == 'Iz':
|
|
||||||
ax = plt.subplot(gs[2, 2])
|
|
||||||
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
|
||||||
ax.set_ylabel(outputtext + ', current [A]')
|
|
||||||
for ax in fig.axes:
|
|
||||||
ax.set_xlim([0, np.amax(time)])
|
|
||||||
ax.grid()
|
|
||||||
|
|
||||||
# Save a PDF/PNG of the figure
|
# Plotting if FFT required
|
||||||
#fig.savefig(os.path.splitext(os.path.abspath(file))[0] + '_rx' + str(rx) + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
if fft:
|
||||||
#fig.savefig(os.path.splitext(os.path.abspath(file))[0] + '_rx' + str(rx) + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
# Calculate magnitude of frequency spectra of waveform
|
||||||
|
power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)
|
||||||
|
freqs = np.fft.fftfreq(power.size, d=dt)
|
||||||
|
|
||||||
plt.show()
|
# Shift powers so that frequency with maximum power is at zero decibels
|
||||||
|
power -= np.amax(power)
|
||||||
|
|
||||||
|
# Set plotting range to -60dB from maximum power
|
||||||
|
pltrange = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
|
||||||
|
# To a maximum frequency
|
||||||
|
#pltrange = np.where(freqs > 2e9)[0][0]
|
||||||
|
pltrange = np.s_[0:pltrange]
|
||||||
|
|
||||||
|
# Plot time history of output component
|
||||||
|
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||||
|
line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext)
|
||||||
|
ax1.set_xlabel('Time [s]')
|
||||||
|
ax1.set_ylabel(outputtext + ' field strength [V/m]')
|
||||||
|
ax1.set_xlim([0, np.amax(time)])
|
||||||
|
ax1.grid()
|
||||||
|
|
||||||
|
# Plot frequency spectra
|
||||||
|
markerline, stemlines, baseline = ax2.stem(freqs[pltrange], power[pltrange], '-.')
|
||||||
|
plt.setp(baseline, 'linewidth', 0)
|
||||||
|
plt.setp(stemlines, 'color', 'r')
|
||||||
|
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||||
|
line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||||
|
ax2.set_xlabel('Frequency [Hz]')
|
||||||
|
ax2.set_ylabel('Power [dB]')
|
||||||
|
ax2.grid()
|
||||||
|
|
||||||
|
# Change colours and labels for magnetic field components or currents
|
||||||
|
if 'H' in outputs[0]:
|
||||||
|
plt.setp(line1, color='g')
|
||||||
|
plt.setp(line2, color='g')
|
||||||
|
plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]')
|
||||||
|
plt.setp(stemlines, 'color', 'g')
|
||||||
|
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||||
|
elif 'I' in outputs[0]:
|
||||||
|
plt.setp(line1, color='b')
|
||||||
|
plt.setp(line2, color='b')
|
||||||
|
plt.setp(ax1, ylabel=outputtext + ' current [A]')
|
||||||
|
plt.setp(stemlines, 'color', 'b')
|
||||||
|
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||||
|
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
# Plotting if no FFT required
|
||||||
|
else:
|
||||||
|
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]', ylabel=outputtext + ' field strength [V/m]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||||
|
line = ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
||||||
|
ax.set_xlim([0, np.amax(time)])
|
||||||
|
#ax.set_ylim([-15, 20])
|
||||||
|
ax.grid()
|
||||||
|
|
||||||
|
if 'H' in output:
|
||||||
|
plt.setp(line, color='g')
|
||||||
|
plt.setp(ax, ylabel=outputtext + ', field strength [A/m]')
|
||||||
|
elif 'I' in output:
|
||||||
|
plt.setp(line, color='b')
|
||||||
|
plt.setp(ax, ylabel=outputtext + ', current [A]')
|
||||||
|
|
||||||
|
# If multiple outputs required, create all nine subplots and populate only the specified ones
|
||||||
|
else:
|
||||||
|
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||||
|
gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3)
|
||||||
|
for output in outputs:
|
||||||
|
|
||||||
|
# Check for polarity of output and if requested output is in file
|
||||||
|
if output[0] == 'm':
|
||||||
|
polarity = -1
|
||||||
|
outputtext = '-' + output[1:]
|
||||||
|
output = output[1:]
|
||||||
|
else:
|
||||||
|
polarity = 1
|
||||||
|
outputtext = output
|
||||||
|
|
||||||
|
# Check if requested output is in file
|
||||||
|
if output not in availableoutputs:
|
||||||
|
raise CmdInputError('Output(s) requested to plot: {}, but available output(s) for receiver {} in the file: {}'.format(', '.join(outputs), rx, ', '.join(availableoutputs)))
|
||||||
|
|
||||||
|
outputdata = f[path + output][:] * polarity
|
||||||
|
|
||||||
|
if output == 'Ex':
|
||||||
|
ax = plt.subplot(gs[0, 0])
|
||||||
|
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
||||||
|
#ax.set_ylim([-15, 20])
|
||||||
|
elif output == 'Ey':
|
||||||
|
ax = plt.subplot(gs[1, 0])
|
||||||
|
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
||||||
|
#ax.set_ylim([-15, 20])
|
||||||
|
elif output == 'Ez':
|
||||||
|
ax = plt.subplot(gs[2, 0])
|
||||||
|
ax.plot(time, outputdata,'r', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [V/m]')
|
||||||
|
#ax.set_ylim([-15, 20])
|
||||||
|
elif output == 'Hx':
|
||||||
|
ax = plt.subplot(gs[0, 1])
|
||||||
|
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
||||||
|
#ax.set_ylim([-0.03, 0.03])
|
||||||
|
elif output == 'Hy':
|
||||||
|
ax = plt.subplot(gs[1, 1])
|
||||||
|
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
||||||
|
#ax.set_ylim([-0.03, 0.03])
|
||||||
|
elif output == 'Hz':
|
||||||
|
ax = plt.subplot(gs[2, 1])
|
||||||
|
ax.plot(time, outputdata,'g', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', field strength [A/m]')
|
||||||
|
#ax.set_ylim([-0.03, 0.03])
|
||||||
|
elif output == 'Ix':
|
||||||
|
ax = plt.subplot(gs[0, 2])
|
||||||
|
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', current [A]')
|
||||||
|
elif output == 'Iy':
|
||||||
|
ax = plt.subplot(gs[1, 2])
|
||||||
|
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', current [A]')
|
||||||
|
elif output == 'Iz':
|
||||||
|
ax = plt.subplot(gs[2, 2])
|
||||||
|
ax.plot(time, outputdata,'b', lw=2, label=outputtext)
|
||||||
|
ax.set_ylabel(outputtext + ', current [A]')
|
||||||
|
for ax in fig.axes:
|
||||||
|
ax.set_xlim([0, np.amax(time)])
|
||||||
|
ax.grid()
|
||||||
|
|
||||||
|
# Save a PDF/PNG of the figure
|
||||||
|
#fig.savefig(os.path.splitext(os.path.abspath(file))[0] + '_rx' + str(rx) + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||||
|
#fig.savefig(os.path.splitext(os.path.abspath(file))[0] + '_rx' + str(rx) + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||||
|
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
|
||||||
|
# Parse command line arguments
|
||||||
|
parser = argparse.ArgumentParser(description='Plots electric and magnetic fields and currents from all receiver points in the given output file. Each receiver point is plotted in a new figure window.', usage='cd gprMax; python -m tools.plot_Ascan outputfile')
|
||||||
|
parser.add_argument('outputfile', help='name of output file including path')
|
||||||
|
parser.add_argument('--outputs', help='outputs to be plotted', default=Rx.availableoutputs, choices='Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, Iz', nargs='+')
|
||||||
|
parser.add_argument('-fft', action='store_true', help='plot FFT (single output must be specified)', default=False)
|
||||||
|
args = parser.parse_args()
|
||||||
|
|
||||||
|
make_plot(args.outputfile, args.outputs, fft=args.fft)
|
||||||
|
@@ -23,52 +23,64 @@ import matplotlib.pyplot as plt
|
|||||||
|
|
||||||
from gprMax.exceptions import CmdInputError
|
from gprMax.exceptions import CmdInputError
|
||||||
|
|
||||||
"""Plots a B-scan image."""
|
|
||||||
|
|
||||||
# Parse command line arguments
|
def make_plot(filename, output):
|
||||||
parser = argparse.ArgumentParser(description='Plots a B-scan image.', usage='cd gprMax; python -m tools.plot_Bscan outputfile output')
|
"""Plots a B-scan image.
|
||||||
parser.add_argument('outputfile', help='name of output file including path')
|
|
||||||
parser.add_argument('output', help='name of output component to be plotted', choices='Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy or Iz')
|
Args:
|
||||||
args = parser.parse_args()
|
filename (string): Filename (including path) of output file.
|
||||||
|
output (string): Field/current component to plot.
|
||||||
|
"""
|
||||||
|
|
||||||
# Open output file and read some attributes
|
# Open output file and read some attributes
|
||||||
f = h5py.File(args.outputfile, 'r')
|
f = h5py.File(filename, 'r')
|
||||||
nrx = f.attrs['nrx']
|
nrx = f.attrs['nrx']
|
||||||
|
|
||||||
# Check there are any receivers
|
# Check there are any receivers
|
||||||
if nrx == 0:
|
if nrx == 0:
|
||||||
raise CmdInputError('No receivers found in {}'.format(args.outputfile))
|
raise CmdInputError('No receivers found in {}'.format(filename))
|
||||||
|
|
||||||
for rx in range(1, nrx + 1):
|
for rx in range(1, nrx + 1):
|
||||||
path = '/rxs/rx' + str(rx) + '/'
|
path = '/rxs/rx' + str(rx) + '/'
|
||||||
availableoutputs = list(f[path].keys())
|
availableoutputs = list(f[path].keys())
|
||||||
|
|
||||||
# Check if requested output is in file
|
# Check if requested output is in file
|
||||||
if args.output not in availableoutputs:
|
if output not in availableoutputs:
|
||||||
raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(args.output, ', '.join(availableoutputs)))
|
raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(output, ', '.join(availableoutputs)))
|
||||||
|
|
||||||
outputdata = f[path + '/' + args.output]
|
outputdata = f[path + '/' + output]
|
||||||
|
|
||||||
# Check that there is more than one A-scan present
|
# Check that there is more than one A-scan present
|
||||||
if outputdata.shape[1] == 1:
|
if outputdata.shape[1] == 1:
|
||||||
raise CmdInputError('{} contains only a single A-scan.'.format(args.outputfile))
|
raise CmdInputError('{} contains only a single A-scan.'.format(filename))
|
||||||
|
|
||||||
# Plot B-scan image
|
# Plot B-scan image
|
||||||
fig = plt.figure(num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
fig = plt.figure(num='rx' + str(rx), figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||||
plt.imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0]*f.attrs['dt'], 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))
|
plt.imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0]*f.attrs['dt'], 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))
|
||||||
plt.xlabel('Trace number')
|
plt.xlabel('Trace number')
|
||||||
plt.ylabel('Time [s]')
|
plt.ylabel('Time [s]')
|
||||||
plt.grid()
|
plt.grid()
|
||||||
cb = plt.colorbar()
|
cb = plt.colorbar()
|
||||||
if 'E' in args.output:
|
if 'E' in args.output:
|
||||||
cb.set_label('Field strength [V/m]')
|
cb.set_label('Field strength [V/m]')
|
||||||
elif 'H' in args.output:
|
elif 'H' in args.output:
|
||||||
cb.set_label('Field strength [A/m]')
|
cb.set_label('Field strength [A/m]')
|
||||||
elif 'I' in args.output:
|
elif 'I' in args.output:
|
||||||
cb.set_label('Current [A]')
|
cb.set_label('Current [A]')
|
||||||
|
|
||||||
# Save a PDF/PNG of the figure
|
# Save a PDF/PNG of the figure
|
||||||
#fig.savefig(os.path.splitext(os.path.abspath(args.outputfile))[0] + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
#fig.savefig(os.path.splitext(os.path.abspath(args.outputfile))[0] + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||||
#fig.savefig(os.path.splitext(os.path.abspath(args.outputfile))[0] + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
#fig.savefig(os.path.splitext(os.path.abspath(args.outputfile))[0] + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||||
|
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
|
||||||
|
# Parse command line arguments
|
||||||
|
parser = argparse.ArgumentParser(description='Plots a B-scan image.', usage='cd gprMax; python -m tools.plot_Bscan outputfile output')
|
||||||
|
parser.add_argument('outputfile', help='name of output file including path')
|
||||||
|
parser.add_argument('output', help='name of output component to be plotted', choices='Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy or Iz')
|
||||||
|
args = parser.parse_args()
|
||||||
|
|
||||||
|
make_plot(args.outputfile, args.output)
|
@@ -25,8 +25,6 @@ from gprMax.exceptions import CmdInputError
|
|||||||
from gprMax.utilities import round_value
|
from gprMax.utilities import round_value
|
||||||
from gprMax.waveforms import Waveform
|
from gprMax.waveforms import Waveform
|
||||||
|
|
||||||
"""Plot built-in waveforms that can be used with sources."""
|
|
||||||
|
|
||||||
|
|
||||||
def check_timewindow(timewindow, dt):
|
def check_timewindow(timewindow, dt):
|
||||||
"""Checks and sets time window and number of iterations.
|
"""Checks and sets time window and number of iterations.
|
||||||
@@ -58,7 +56,7 @@ def check_timewindow(timewindow, dt):
|
|||||||
return timewindow, iterations
|
return timewindow, iterations
|
||||||
|
|
||||||
|
|
||||||
def plot_waveform(w, timewindow, dt, iterations, fft=False):
|
def make_plot(w, timewindow, dt, iterations, fft=False):
|
||||||
"""Plots waveform and prints useful information about its properties.
|
"""Plots waveform and prints useful information about its properties.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
@@ -66,7 +64,7 @@ def plot_waveform(w, timewindow, dt, iterations, fft=False):
|
|||||||
timewindow (float): Time window.
|
timewindow (float): Time window.
|
||||||
dt (float): Time discretisation.
|
dt (float): Time discretisation.
|
||||||
iterations (int): Number of iterations.
|
iterations (int): Number of iterations.
|
||||||
fft (boolean): Plot FFT of waveform.
|
fft (boolean): Plot FFT switch.
|
||||||
"""
|
"""
|
||||||
|
|
||||||
time = np.linspace(0, 1, iterations)
|
time = np.linspace(0, 1, iterations)
|
||||||
@@ -118,7 +116,7 @@ def plot_waveform(w, timewindow, dt, iterations, fft=False):
|
|||||||
plt.setp(baseline, 'linewidth', 0)
|
plt.setp(baseline, 'linewidth', 0)
|
||||||
plt.setp(stemlines, 'color', 'r')
|
plt.setp(stemlines, 'color', 'r')
|
||||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||||
ax2.plot(freqs[pltrange]/1e9, power[pltrange], 'r', lw=2)
|
ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2)
|
||||||
ax2.set_xlabel('Frequency [Hz]')
|
ax2.set_xlabel('Frequency [Hz]')
|
||||||
ax2.set_ylabel('Power [dB]')
|
ax2.set_ylabel('Power [dB]')
|
||||||
|
|
||||||
@@ -148,7 +146,7 @@ if __name__ == "__main__":
|
|||||||
parser.add_argument('freq', type=float, help='centre frequency of waveform')
|
parser.add_argument('freq', type=float, help='centre frequency of waveform')
|
||||||
parser.add_argument('timewindow', help='time window to view waveform')
|
parser.add_argument('timewindow', help='time window to view waveform')
|
||||||
parser.add_argument('dt', type=float, help='time step to view waveform')
|
parser.add_argument('dt', type=float, help='time step to view waveform')
|
||||||
parser.add_argument('-fft', action='store_true', default=False, help='plot FFT of waveform')
|
parser.add_argument('-fft', action='store_true', help='plot FFT of waveform', default=False)
|
||||||
args = parser.parse_args()
|
args = parser.parse_args()
|
||||||
|
|
||||||
# Check waveform parameters
|
# Check waveform parameters
|
||||||
@@ -164,7 +162,7 @@ if __name__ == "__main__":
|
|||||||
w.freq = args.freq
|
w.freq = args.freq
|
||||||
|
|
||||||
timewindow, iterations = check_timewindow(args.timewindow, args.dt)
|
timewindow, iterations = check_timewindow(args.timewindow, args.dt)
|
||||||
plot_waveform(w, timewindow, args.dt, iterations, args.fft)
|
make_plot(w, timewindow, args.dt, iterations, args.fft)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
在新工单中引用
屏蔽一个用户