Obspy: how to read the files by one and one ?

[Earthscience]

OSX

import numpy as np
import matplotlib.pyplot as plt
import obspy
import os
import sys
from obspy.core import read
from obspy.signal.tf_misfit import cwt
import pylab, pysac
import time
import glob


for hour in range(0,24):
    for minute in range(0,6):

          if minute < 6: #10~60 minute seqeunce


                    if hour <= 9:
                              st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,%s%s:%s0:00.SAC" %(0, hour, minute))
                              print st
                              tr = st[0]
                              continue


                        elif hour >= 10:
                            st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,%s:%s0:00.SAC" %(hour, minute))



                            if hour == 24:
                                    break

                            print st

and add spectrogram code.

i want to read the file sequentially, want to be spectrogram code is executed.
but this code not working.

This code reads only the last file.

how to do problem?

[krischer]

I don't really understand what you want to do to be honest but maybe this helps: you don't need the conditional - you can just always do this:

st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,%02i:%i0:00.SAC" %(hour, minute))

because

In [1]: "%02i" % 2
Out[1]: '02'

In [2]: "%02i" % 10
Out[2]: '10'

[trichter]

Without the full code it is difficult to guess what is wrong. You are reading in file by file, maybe you forget to add the streams together? You can also read in all files at once into a single stream with e.g.

st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,*:*0:00.SAC")

or maybe just

st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.SAC")

[claudiodsf]

I took the liberty of reformatting your code for better readability:

import numpy as np
import matplotlib.pyplot as plt
import obspy
import os
import sys
from obspy.core import read
from obspy.signal.tf_misfit import cwt
import pylab
import pysac
import time
import glob

for hour in range(0, 24):
    for minute in range(0, 6):

        if minute < 6:  # 10~60 minute seqeunce

            if hour <= 9:

                st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,*:*0:00.SAC")
                print st
                continue

            elif hour >= 10:
                st = read(
                    "/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,%s:%s0:00.SAC" %
                    (hour, minute))

                if hour == 24:
                    break

                print st


#tr = read("/Users/EE/Desktop/03/05/TP.TPJ02..HHZ.D.2015,064,%s0:%s0:00.SAC" %(hour, tt), format="SAC")[0]
tr = st[0]


tr.detrend('linear')
tr.detrend('demean')

# bandpass 1-5hz
#tr.filter("bandpass", freqmin=1, freqmax=5)

# file name edit
ee = str(tr)

npts = tr.stats.npts
dt = tr.stats.delta
t = np.linspace(0, dt * npts, npts)
f_min = 1
f_max = 20

scalogram = cwt(tr.data, dt, 8, f_min, f_max)

fig = plt.figure()
ax1 = fig.add_axes([0.09, 0.1, 0.80, 0.60])
ax2 = fig.add_axes([0.09, 0.75, 0.80, 0.2])
ax3 = fig.add_axes([0.91, 0.1, 0.03, 0.6])

img = ax1.imshow(np.abs(scalogram)[-1::-1],
                 extent=[t[0], t[-1], f_min, f_max],
                 aspect='auto',
                 interpolation="nearest")

ax1.set_xlabel("Time after %s [s]" % tr.stats.starttime, fontsize=18)
ax1.set_ylabel("Frequency [Hz]", fontsize=15)
ax1.set_yscale('linear', fontsize=5)
ax2.plot(t, tr.data, 'k')

pylab.xlim([0, 400])


fig.colorbar(img, cax=ax3)

# size change
fig.set_size_inches(16, 8)


# save figure directory, dpi
#plt.savefig('11.png', dpi=400, bbox_inches='tight')
fig.savefig('/Users/EE/Desktop/Spectrogram/03/06/' +
            ee[16:35] + 'Z.png', dpi=200)


# show figure
# plt.show()

[claudiodsf]

One problem I see is that you overwrite st at every loop.

To avoid this you can declare an empty Stream object:

from obspy import Stream
st = Stream()

and use the += operator to append to it:

st += read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,*:*0:00.SAC")

st += read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,%s:%s0:00.SAC" % (hour, minute))

[claudiodsf]

More easily, you can read all the files without any loop, by using something like this:

st = read("/Users/EE/Desktop/03/05/TP.TPJ02.*.HHZ*,[0-2][0-9]:[0-6]0:00.SAC")

Each bracket (e.g., [0-2]) is what we call a "regular expression". It gives an interval of variability for that particular character. Therefore [0-2][0-9]:[0-6] means that:

• [0-2] the hour first digit can vary between 0 and 2
• [0-9] the hour second digit can vary between 0 and 9
• [0-6] the minute first digit can vary between 0 and 6

[claudiodsf]

You can loop over traces in stream with for tr in st:
Therefore your code will become:

for tr in st:

    tr.detrend('linear')
    tr.detrend('demean')

    # bandpass 1-5hz
    #tr.filter("bandpass", freqmin=1, freqmax=5)

    # file name edit
    ee = str(tr)

    npts = tr.stats.npts
    dt = tr.stats.delta
    t = np.linspace(0, dt * npts, npts)
    f_min = 1
    f_max = 20

    scalogram = cwt(tr.data, dt, 8, f_min, f_max)

    fig = plt.figure()
    ax1 = fig.add_axes([0.09, 0.1, 0.80, 0.60])
    ax2 = fig.add_axes([0.09, 0.75, 0.80, 0.2])
    ax3 = fig.add_axes([0.91, 0.1, 0.03, 0.6])

    img = ax1.imshow(np.abs(scalogram)[-1::-1],
                     extent=[t[0], t[-1], f_min, f_max],
                     aspect='auto',
                     interpolation="nearest")

    ax1.set_xlabel("Time after %s [s]" % tr.stats.starttime, fontsize=18)
    ax1.set_ylabel("Frequency [Hz]", fontsize=15)
    ax1.set_yscale('linear', fontsize=5)
    ax2.plot(t, tr.data, 'k')

    pylab.xlim([0, 400])


    fig.colorbar(img, cax=ax3)

    # size change
    fig.set_size_inches(16, 8)


    # save figure directory, dpi
    #plt.savefig('11.png', dpi=400, bbox_inches='tight')
    fig.savefig('/Users/EE/Desktop/Spectrogram/03/06/' +
                ee[16:35] + 'Z.png', dpi=200)

[Earthscience]

thanks! very very thank you!!..

modified code is operated successfully.

[QuLogic]

These sorts of questions are probably better off on the mailing list, but good you got it working.