Python core.Trace类代码示例

def export_sac(db, filename, pair, components, filterid, corr, ncorr=0, sac_format=None, maxlag=None, cc_sampling_rate=None):
    if sac_format is None:
        sac_format = get_config(db, "sac_format")
    if maxlag is None:
        maxlag = float(get_config(db, "maxlag"))
    if cc_sampling_rate is None:
        cc_sampling_rate = float(get_config(db, "cc_sampling_rate"))
    try:
        os.makedirs(os.path.split(filename)[0])
    except:
        pass
    filename += ".SAC"
    mytrace = Trace(data=corr)
    mytrace.stats['station'] = pair
    mytrace.stats['sampling_rate'] = cc_sampling_rate

    st = Stream(traces=[mytrace, ])
    st.write(filename, format='SAC')
    tr = SacIO(filename)
    if sac_format == "doublets":
        tr.SetHvalue('A', 120)
    else:
        tr.SetHvalue('B', -maxlag)
        tr.SetHvalue('DEPMIN', np.min(corr))
        tr.SetHvalue('DEPMAX', np.max(corr))
        tr.SetHvalue('DEPMEN', np.mean(corr))
        tr.SetHvalue('SCALE', 1)
        tr.SetHvalue('NPTS', len(corr))
    tr.WriteSacBinary(filename)
    del st, tr
    return

def export_sac(db, filename, pair, components, filterid, corr, ncorr=0,
               sac_format=None, maxlag=None, cc_sampling_rate=None):
    if sac_format is None:
        sac_format = get_config(db, "sac_format")
    if maxlag is None:
        maxlag = float(get_config(db, "maxlag"))
    if cc_sampling_rate is None:
        cc_sampling_rate = float(get_config(db, "cc_sampling_rate"))
    try:
        os.makedirs(os.path.split(filename)[0])
    except:
        pass
    filename += ".SAC"
    mytrace = Trace(data=corr)
    mytrace.stats['station'] = pair
    mytrace.stats['sampling_rate'] = cc_sampling_rate
    if maxlag:
        mytrace.stats.starttime = -maxlag
    mytrace.stats.sac = AttribDict()
    mytrace.stats.sac.depmin = np.min(corr)
    mytrace.stats.sac.depmax = np.max(corr)
    mytrace.stats.sac.depmen = np.mean(corr)
    mytrace.stats.sac.scale = 1
    mytrace.stats.sac.npts = len(corr)

    st = Stream(traces=[mytrace, ])
    st.write(filename, format='SAC')
    del st
    return

 def test_len(self):
     """
     Tests the __len__ and count methods of the Trace class.
     """
     trace = Trace(data=np.arange(1000))
     self.assertEquals(len(trace), 1000)
     self.assertEquals(trace.count(), 1000)

 def setUp(self):
     # directory where the test files are located
     self.path = os.path.join(os.path.dirname(__file__), 'data')
     self.filename_css = os.path.join(self.path, 'test_css.wfdisc')
     self.filename_nnsa = os.path.join(self.path, 'test_nnsa.wfdisc')
     # set up stream for validation
     header = {}
     header['station'] = 'TEST'
     header['starttime'] = UTCDateTime(1296474900.0)
     header['sampling_rate'] = 80.0
     header['calib'] = 1.0
     header['calper'] = 1.0
     header['_format'] = 'CSS'
     filename = os.path.join(self.path, '201101311155.10.ascii.gz')
     with gzip.open(filename, 'rb') as fp:
         data = np.loadtxt(fp, dtype=np.int_)
     # traces in the test files are sorted ZEN
     st = Stream()
     for x, cha in zip(data.reshape((3, 4800)), ('HHZ', 'HHE', 'HHN')):
         # big-endian copy
         tr = Trace(x, header.copy())
         tr.stats.station += 'be'
         tr.stats.channel = cha
         st += tr
         # little-endian copy
         tr = Trace(x, header.copy())
         tr.stats.station += 'le'
         tr.stats.channel = cha
         st += tr
     self.st_result_css = st.copy()
     for tr in st:
         tr.stats['_format'] = "NNSA_KB_CORE"
     self.st_result_nnsa = st

 def test_percent_in_str(self):
     """
     Tests if __str__ method is working with percent sign (%).
     """
     tr = Trace()
     tr.stats.station = '%t3u'
     self.assertTrue(tr.__str__().startswith(".%t3u.. | 1970"))

 def test_issue193(self):
     """
     Test for issue #193: if non-contiguous array is written correctly.
     """
     warnings.filterwarnings("ignore", "Detected non contiguous data")
     # test all plugins with both read and write method
     formats_write = \
         set(_getEntryPoints('obspy.plugin.waveform', 'writeFormat'))
     formats_read = \
         set(_getEntryPoints('obspy.plugin.waveform', 'readFormat'))
     formats = set.intersection(formats_write, formats_read)
     # mseed will raise exception for int64 data, thus use int32 only
     data = np.arange(10, dtype='int32')
     # make array non-contiguous
     data = data[::2]
     tr = Trace(data=data)
     for format in formats:
         # XXX: skip SEGY and SU formats for now as they need some special
         # headers.
         if format in ['SEGY', 'SU', 'SEG2']:
             continue
         tempfile = NamedTemporaryFile().name
         tr.write(tempfile, format)
         if format == "Q":
             tempfile = tempfile + ".QHD"
         tr_test = read(tempfile, format)[0]
         # clean up
         os.remove(tempfile)
         if format == 'Q':
             os.remove(tempfile[:-4] + '.QBN')
             os.remove(tempfile[:-4])
         np.testing.assert_array_equal(tr.data, tr_test.data)

    def test_SacInstCorrection(self):
        # SAC recommends to taper the transfer function if a pure
        # deconvolution is done instead of simulating a different
        # instrument. This test checks the difference between the
        # result from removing the instrument response using SAC or
        # ObsPy. Visual inspection shows that the traces are pretty
        # much identical but differences remain (rms ~ 0.042). Haven't
        # found the cause for those, yet. One possible reason is the
        # floating point arithmetic of SAC vs. the double precision
        # arithmetic of Python. However differences still seem to be
        # too big for that.
        pzf = os.path.join(self.path, 'SAC_PZs_KARC_BHZ')
        sacf = os.path.join(self.path, 'KARC.LHZ.SAC.asc.gz')
        testsacf = os.path.join(self.path, 'KARC_corrected.sac.asc.gz')
        plow = 160.
        phigh = 4.
        fl1 = 1.0 / (plow + 0.0625 * plow)
        fl2 = 1.0 / plow
        fl3 = 1.0 / phigh
        fl4 = 1.0 / (phigh - 0.25 * phigh)
        #Uncomment the following to run the sac-commands
        #that created the testing file
        #if 1:
        #    import subprocess as sp
        #    p = sp.Popen('sac',shell=True,stdin=sp.PIPE)
        #    cd1 = p.stdin
        #    print >>cd1, "r %s"%sacf
        #    print >>cd1, "rmean"
        #    print >>cd1, "rtrend"
        #    print >>cd1, "taper type cosine width 0.03"
        #    print >>cd1, "transfer from polezero subtype %s to none \
        #    freqlimits %f %f %f %f" % (pzf, fl1, fl2, fl3, fl4)
        #    print >>cd1, "w over ./data/KARC_corrected.sac"
        #    print >>cd1, "quit"
        #    cd1.close()
        #    p.wait()

        stats = {'network': 'KA', 'delta': 0.99999988079072466,
                 'station': 'KARC', 'location': 'S1',
                 'starttime': UTCDateTime(2001, 2, 13, 0, 0, 0, 993700),
                 'calib': 1.00868e+09, 'channel': 'BHZ'}
        tr = Trace(np.loadtxt(sacf), stats)

        attach_paz(tr, pzf, tovel=False)
        tr.data = seisSim(tr.data, tr.stats.sampling_rate,
                          paz_remove=tr.stats.paz, remove_sensitivity=False,
                          pre_filt=(fl1, fl2, fl3, fl4))

        data = np.loadtxt(testsacf)

        # import matplotlib.pyplot as plt
        # plt.plot(tr.data)
        # plt.plot(data)
        # plt.show()
        rms = np.sqrt(np.sum((tr.data - data) ** 2) / \
                      np.sum(tr.data ** 2))
        self.assertTrue(rms < 0.0421)

 def test_integrate(self):
     """
     Test integration method of trace
     """
     data = np.ones(101) * 0.01
     tr = Trace(data=data)
     tr.stats.delta = 0.1
     tr.integrate(type='cumtrapz')
     np.testing.assert_almost_equal(tr.data[-1], 0.1)

 def test_taper(self):
     """
     Test taper method of trace
     """
     data = np.ones(10)
     tr = Trace(data=data)
     tr.taper()
     for i in range(len(data)):
         self.assertTrue(tr.data[i] <= 1.)
         self.assertTrue(tr.data[i] >= 0.)

 def test_differentiate(self):
     """
     Test differentiation method of trace
     """
     t = np.linspace(0., 1., 11)
     data = 0.1 * t + 1.
     tr = Trace(data=data)
     tr.stats.delta = 0.1
     tr.differentiate(type='gradient')
     np.testing.assert_array_almost_equal(tr.data, np.ones(11) * 0.1)

 def test_taper(self):
     """
     Test taper method of trace
     """
     data = np.ones(10)
     tr = Trace(data=data)
     tr.taper()
     for i in range(len(data)):
         self.assertLessEqual(tr.data[i], 1.)
         self.assertGreaterEqual(tr.data[i], 0.)

def semblance(st, s, baz, winlen):
    """
    Returns the semblance for a seismic array, for a beam of given slowness and backazimuth.

    Parameters
    ----------
    st : ObsPy Stream object
        Stream of SAC format seismograms for the seismic array, length K = no. of stations in array
    s  : float
        Magnitude of slowness vector, in s / km
    baz : float
        Backazimuth of slowness vector, (i.e. angle from North back to epicentre of event)
    winlen : int
        Length of Hann window over which to calculate the semblance.

    Returns
    -------
    semblance : NumPy array
        The semblance at the given slowness and backazimuth, as a time series.

    """

    # Check that each channel has the same number of samples, otherwise we can't construct the beam properly
    assert len(set([len(tr) for tr in st])) == 1, "Traces in stream have different lengths, cannot stack."

    nsta = len(st)

    stack = linear_stack(st, s, baz)

    # Taper the linear stack
    stack_trace = Trace(stack)
    stack_trace.taper(type="cosine", max_percentage=0.05)
    stack = stack_trace.data

    shifts = get_shifts(st, s, baz)

    # Smooth data with sliding Hann window (i.e. convolution of signal and window function)
    window = np.hanning(winlen)

    shifted_st = st.copy()

    for i, tr in enumerate(shifted_st):
        tr.data = np.roll(tr.data, shifts[i])  # Shift data in each trace by its offset
        tr.taper(type="cosine", max_percentage=0.05)  # Taper it

    # Calculate the power in the beam
    beampower = np.convolve(stack ** 2, window, mode="same")

    # Calculate the summed power of each trace
    tracepower = np.convolve(np.sum([tr.data ** 2 for tr in shifted_st], axis=0), window, mode="same")

    # Calculate semblance
    semblance = nsta * beampower / tracepower

    return semblance

 def test_slice(self):
     """
     Tests the slicing of trace objects.
     """
     tr = Trace(data=np.arange(10, dtype='int32'))
     mempos = tr.data.ctypes.data
     t = tr.stats.starttime
     tr1 = tr.slice(t + 2, t + 8)
     tr1.data[0] = 10
     self.assertEqual(tr.data[2], 10)
     self.assertEqual(tr.data.ctypes.data, mempos)
     self.assertEqual(tr.data[2:9].ctypes.data, tr1.data.ctypes.data)
     self.assertEqual(tr1.data.ctypes.data - 8, mempos)

 def test_writeSACXYWithMinimumStats(self):
     """
     Write SACXY with minimal stats header, no inhereted from SAC file
     """
     tr = Trace()
     tr.stats.delta = 0.01
     tr.data = np.arange(0, 3000)
     sac_file = NamedTemporaryFile().name
     tr.write(sac_file, "SACXY")
     st = read(sac_file)
     os.remove(sac_file)
     self.assertEquals(st[0].stats.delta, 0.01)
     self.assertEquals(st[0].stats.sampling_rate, 100.0)

    def test_trimAllDoesNotChangeDtype(self):
        """
        If a Trace is completely trimmed, e.g. no data samples are remaining,
        the dtype should remain unchanged.

        A trace with no data samples is not really senseful but the dtype
        should not be changed anyways.
        """
        # Choose non native dtype.
        tr = Trace(np.arange(100, dtype='int16'))
        tr.trim(UTCDateTime(10000), UTCDateTime(20000))
        # Assert the result.
        self.assertEqual(len(tr.data), 0)
        self.assertEqual(tr.data.dtype, 'int16')