diff --git a/utils/plot_ASDF_data_info.py b/utils/plot_ASDF_data_info.py new file mode 100755 index 000000000..48e4a03cd --- /dev/null +++ b/utils/plot_ASDF_data_info.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python +# +# see: http://seismicdata.github.io/pyasdf/tutorial.html +# +from __future__ import print_function +import sys +import pyasdf + +def plot_ASDF_data_info(file,show_plot): + print("") + print("reading file: ",filename) + print("") + + ds = pyasdf.ASDFDataSet(filename) + + # file info + print("file info...") + print(ds) + print("") + + # validate + print("validate...") + ds.validate() + print("") + + # event info + print("event info...") + print(" number of events: ",len(ds.events)) + print("") + type(ds.events) + print(ds.events) + for event in ds.events: + print(event) + print("") + + # station info + print("") + print("station list...") + print(" number of stations: ",len(ds.waveforms)) + print("") + print(ds.waveforms.list()) + print("") + + for station in ds.waveforms.list(): + print("") + print("station:") + print("") + print(station) + + # how to get station name and associated waveform? + # waveforms + sta = ds.waveforms[station] + type(sta) + print(sta.synthetic) + print("") + print("waveform list:",sta.list()) + print("waveform tags:",sta.get_waveform_tags()) + print("") + + # plotting + if show_plot: + print("plotting stream...") + st = sta.synthetic + st.plot() + + +def usage(): + print("usage: ./read_asdf_data_info.py filename[e.g. synthetics.h5] [show]") + print(" with") + print(" filename - e.g. synthetics.h5") + print(" show - (optional) plot waveforms") + sys.exit(1) + + +if __name__ == "__main__": + # gets arguments + do_plot_waveforms = False + + if len(sys.argv) == 2: + filename = sys.argv[1] + elif len(sys.argv) == 3: + filename = sys.argv[1] + if sys.argv[2] == "show": do_plot_waveforms = True + else: + usage() + sys.exit(1) + + plot_ASDF_data_info(filename,do_plot_waveforms) diff --git a/utils/plot_SU_file_with_obspy.py b/utils/plot_SU_file_with_obspy.py new file mode 100755 index 000000000..da4df407d --- /dev/null +++ b/utils/plot_SU_file_with_obspy.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python +# +# plot_SU_file_with_obspy.py +# +from __future__ import print_function +import os,sys + +## matplotlib +#import matplotlib as mpl +#print("matplotlib version: ",mpl.__version__) + +import matplotlib.pyplot as plt + +plt.style.use('ggplot') +plt.rcParams['figure.figsize'] = 20, 5 +plt.rcParams['lines.linewidth'] = 0.5 + +## numpy +import numpy as np +print("numpy version: ",np.__version__) + +## do not show scipy warnings +#import warnings +#warnings.filterwarnings('ignore') + +## obspy +import obspy +print("obspy version: ",obspy.__version__) + + +## Seismic Unix +# +# needed on mcmillan +#suoldtonew < $file | suxwigb #perc=99 & + +# reads header info: +# ns = number of samples +# dt = time step +# sx,sy = source coordinates +# gx,gy = receiver coordinates +#sugethw <$sufile key=dt,ns,sx,sy,gx,gy verbose=1 +# +# removes headers +#sustrip <$sufile > $sufile.raw +# +# wiggle plot +#xwigb <$sufile n1=$NSTEP n2=$ntraces & +# +# image plot +#ximage <$sufile n1=$NSTEP n2=$ntraces & +# +# ps image +#suoldtonew <$sufile | supsimage label1="Time (s)" label2="Offset (m)" title="$sufile" perc=99 labelsize=22 verbose=1 > tmp.eps +# +# convert to png +#ps2pdf -r300 -sDEVICE=png16m -sOutputFile=$sufile.png tmp.eps + + + +def plot_SU_file(file,show_plot): + print("") + print("getting station information...") + print("") + + st = obspy.read(file) + + print("traces: ") + print(st) + + print("") + print("number of traces: ",len(st)) + print("") + + npts = 0 + duration = 0.0 + + # trace statistics + if len(st) > 0: + tr = st[0] + #print(tr.stats) + npts = tr.stats.npts # number of samples + dt = tr.stats.delta # time step + duration = npts * dt + + freq_Ny = 1.0/(2.0 * dt) # Nyquist frequency + freq_min = 1./duration # minimal possible frequency for length of trace + + print("trace info:") + print(" number of samples = ",npts) + print(" time step = ",dt,"(s)") + print(" duration = ",duration,"(s)") + print(" Nyquist frequency = ",freq_Ny) + print(" minimum frequency = ",freq_min) + print("") + else: + print("no trace information, exiting...") + sys.exit(1) + + # plotting + if show_plot: + # plots all traces + st.plot() + + # time axis for plotting (in s) + t_d = np.linspace(0,duration,npts) + + # vertical trace + for i,tr in enumerate(st): + print("Trace "+str(i+1)+":") + print(tr) + #print(tr.stats) + # plotting + plt.title("Trace " + str(i+1)) + # vertical component + plt.plot(t_d, tr.data, color='black', linewidth=1.0,label="raw data trace "+str(i+1)) + plt.legend() + plt.xlabel("Time (s)") + plt.ylabel("Counts") + plt.show() + + +def usage(): + print("usage: ./plot_SU_file_with_obspy.py su_file [show]") + print(" with") + print(" su_file - seismogram file e.g. OUTPUT_FILES/0_dz_SU") + print(" show - (optional) plot waveforms") + sys.exit(1) + + +if __name__ == "__main__": + # gets arguments + do_plot_waveforms = False + if len(sys.argv) == 2: + file = sys.argv[1] + elif len(sys.argv) == 3: + file = sys.argv[1] + if sys.argv[2] == "show": do_plot_waveforms = True + else: + usage() + sys.exit(1) + + plot_SU_file(file,do_plot_waveforms) + + + +