Python core.Stream类代码示例

 def test_writeIntegersViaObsPy(self):
     """
     Write file test via L{obspy.Trace}.
     """
     tempfile = NamedTemporaryFile().name
     npts = 1000
     # data cloud of integers - float won't work!
     np.random.seed(815)  # make test reproducable
     data = np.random.randint(-1000, 1000, npts).astype('int32')
     stats = {'network': 'BW', 'station': 'TEST', 'location': '',
              'channel': 'EHE', 'npts': npts, 'sampling_rate': 200.0}
     start = UTCDateTime(2000, 1, 1)
     stats['starttime'] = start
     tr = Trace(data=data, header=stats)
     st = Stream([tr])
     st.verify()
     # write
     st.write(tempfile, format="GSE2")
     # read again
     stream = read(tempfile)
     os.remove(tempfile)
     stream.verify()
     np.testing.assert_equal(data, stream[0].data)
     # test default attributes
     self.assertEqual('CM6', stream[0].stats.gse2.datatype)
     self.assertEqual(-1, stream[0].stats.gse2.vang)
     self.assertEqual(1.0, stream[0].stats.gse2.calper)
     self.assertEqual(1.0, stream[0].stats.calib)

def get_waveforms(session, wfdisc, station=None, channel=None, starttime=None, 
        endtime=None, wfids=None):
    """
    Get waveforms.

    Parameters
    ----------
    session : sqlalchemy.orm.Session instance
        Must be bound.
    wfdisc : mapped Wfdisc table class
    station, channel : str, optional
        Desired station, channel code strings
    starttimes, endtimes : float, optional
        Epoch start times, end times.  
        Traces will be cut to these times.
    wfids : iterable of int, optional
        Wfdisc wfids.  Obviates the above arguments and just returns full Wfdisc
        row waveforms.

    Returns
    -------
    obspy.Stream
        Traces are merged and cut to requested times.

    """
    #TODO: add evids= option?, use with stawin= option in .execute method?
    #TODO: implement get_arrivals if arrivals=True
    Wfdisc = wfdisc

    st = Stream()
    if not wfids:
        t1 = float(starttime)
        t2 = float(endtime)
        sta = station
        chan = channel

        t1_utc = UTCDateTime(float(t1))
        t2_utc = UTCDateTime(float(t2))

        wfs = get_wfdisc_rows( session,Wfdisc, sta, chan, t1, t2)

        #TODO: do arrival stuff here
        for wf in wfs:
            try:
                tr = wfdisc2trace(wf)
                tr.trim(t1_utc, t2_utc)
                st.append(tr)
            except AttributeError:
                #tr is None b/c data couldn't be read
                pass
    else:
        wfs = get_wfdisc_rows( session,Wfdisc, wfids=wfids)
        for wf in wfs:
            try:
                tr = wfdisc2trace(wf)
                st.append(tr)
            except AttributeError:
                pass

    return st

def axisem2mseed_all(path):
    """
    change .dat files into MSEED format
    """
    global test_param

    t = UTCDateTime(0)
    traces = []

    for file in glob.iglob(os.path.join(path, '*.dat')):
        stationID = file.split('/')[-1].split('_')[0]
        chan = file.split('/')[-1].split('_')[-1].split('.')[0]
        dat = np.loadtxt(file)
        npts = len(dat[:,0])
        stats = {'network': 'SG',
                 'station': stationID,
                 'location': '',
                 'channel': chan,
                 'npts': npts,
                 'sampling_rate': (npts - 1.)/(dat[-1,0] - dat[0,0]),
                 'starttime': t + dat[0,0],
                 'mseed' : {'dataquality': 'D'}}
        traces.append(Trace(data=dat[:,1], header=stats))

    st = Stream(traces)
    st.sort()
    fname =  os.path.join(path, 'seismograms.mseed')
    print fname
    st.write(fname, format='MSEED')

 def test_writeAndReadDifferentRecordLengths(self):
     """
     Tests Mini-SEED writing and record lengths.
     """
     # libmseed instance.
     npts = 6000
     np.random.seed(815)  # make test reproducable
     data = np.random.randint(-1000, 1000, npts).astype('int32')
     st = Stream([Trace(data=data)])
     record_lengths = [256, 512, 1024, 2048, 4096, 8192]
     # Loop over some record lengths.
     for rec_len in record_lengths:
         # Write it.
         tempfile = NamedTemporaryFile().name
         st.write(tempfile, format="MSEED", reclen=rec_len)
         # Get additional header info
         info = util.getRecordInformation(tempfile)
         # Test reading the two files.
         temp_st = read(tempfile)
         np.testing.assert_array_equal(data, temp_st[0].data)
         del temp_st
         os.remove(tempfile)
         # Check record length.
         self.assertEqual(info['record_length'], rec_len)
         # Check if filesize is a multiple of the record length.
         self.assertEqual(info['filesize'] % rec_len, 0)

 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 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 test_mergePreviews2(self):
     """
     Test case for issue #84.
     """
     # Note: explicitly creating np.ones instead of np.empty in order to
     # prevent NumPy warnings related to max function
     tr1 = Trace(data=np.ones(2880))
     tr1.stats.starttime = UTCDateTime("2010-01-01T00:00:00.670000Z")
     tr1.stats.delta = 30.0
     tr1.stats.preview = True
     tr1.verify()
     tr2 = Trace(data=np.ones(2881))
     tr2.stats.starttime = UTCDateTime("2010-01-01T23:59:30.670000Z")
     tr2.stats.delta = 30.0
     tr2.stats.preview = True
     tr2.verify()
     st1 = Stream([tr1, tr2])
     st1.verify()
     # merge
     st2 = mergePreviews(st1)
     st2.verify()
     # check
     self.assertTrue(st2[0].stats.preview)
     self.assertEqual(st2[0].stats.starttime, tr1.stats.starttime)
     self.assertEqual(st2[0].stats.endtime, tr2.stats.endtime)
     self.assertEqual(st2[0].stats.npts, 5760)
     self.assertEqual(len(st2[0]), 5760)

def merge_single(nch,dstart,dend):
  '''Merges traces of one channel to larger traces. Used for cross-correlation'''

  # here you load all the functions you need to use
  from obspy.seg2.seg2 import readSEG2
  from obspy.core import Stream


  dataDir2 = "/import/neptun-radler/STEINACH_feb/"
  dataDir = "/import/three-data/hadzii/STEINACH/STEINACH_longtime/"
  outdir = "/home/jsalvermoser/Desktop/Processing/out_merged"

  tr = []

  for k in range(dstart, dend, 1):
    fname = '%d' %(k)
    fileName = fname + ".dat" 
    st = readSEG2(dataDir + fileName)
    #st.detrend('linear')
    tr.append(st[nch-1])


  
  new_stream = Stream(traces=tr)
  new_stream.merge(method=1, fill_value='interpolate')

  start = new_stream[0].stats.starttime
  end = new_stream[0].stats.endtime

  timeframe = str(nch)+ "_" + str(start.year) +'.'+ str(start.julday) +'.'+ str(start.hour) +'.'+ str(start.minute) +'.'+ str(start.second) \
      +'-'+ str(end.year) +'.'+ str(end.julday) +'.'+ str(end.hour) +'.'+ str(end.minute) +'.'+ str(end.second)

  new_stream.write(outdir + timeframe + ".mseed", format="MSEED")

  return new_stream[0]

def ascii2sac(sismograma):
    Wav, Head = [], []
    r=open(sismograma, 'rU')
    
    counter = 0
    for linea in r:
        counter+=1
        if counter <= 5:
            #               print linea
            Head.append(linea.strip())
        else:
            Wav.append(linea)
#print Head

    date = Head[2][0:10]
    hour = Head[2][11:19]
    Fs = float(Head[3][0:8])
    counts = float(Head[4][0:4])
    station = Head[1]
    #si la estacion es triaxial:
    channel = Head[1][3]
    esta = station[0:3]
    
    print date, hour, str(Fs), str(counts), station
    
    data = np.array(Wav, dtype=np.float32)
    #    print data
    
    if channel == 'Z':
        stats = {'network': 'OP', 'station': esta , 'location': '',
        'channel': 'BHZ', 'npts': len(data), 'sampling_rate': Fs,
        'mseed': {'dataquality': 'D'}}
    elif channel == 'N':
        stats = {'network': 'OP', 'station': esta , 'location': '',
        'channel': 'BHN', 'npts': len(data), 'sampling_rate': Fs,
        'mseed': {'dataquality': 'D'}}
    elif channel == 'E':
        stats = {'network': 'OP', 'station': esta , 'location': '',
        'channel': 'BHE', 'npts': len(data), 'sampling_rate': Fs,
        'mseed': {'dataquality': 'D'}}
    else:
        stats = {'network': 'OP', 'station': station , 'location': '',
        'channel': 'SHZ', 'npts': len(data), 'sampling_rate': Fs,
        'mseed': {'dataquality': 'D'}}

    Date = date+hour
    Date = datetime.datetime.strptime(Date, '%Y/%m/%d%H:%M:%S')
    starttime = UTCDateTime(Date)
    stats['starttime'] = starttime

    print stats
    
    st = Stream([Trace(data=data, header=stats)])
    st.write(r.name[0:-4]+'.sac', format='SAC', encoding=4)
    #    st1 = read(r.name[0:-4]+'.mseed')
    #    print st1
    #    st1.plot(color='r')
    return('traza convertida')

#ascii2sac(sismo)

 def test_Header(self):
     """
     Tests whether the header is correctly written and read.
     """
     tempfile = NamedTemporaryFile().name
     np.random.seed(815)  # make test reproducable
     data = np.random.randint(-1000, 1000, 50).astype('int32')
     stats = {'network': 'BW', 'station': 'TEST', 'location': 'A',
              'channel': 'EHE', 'npts': len(data), 'sampling_rate': 200.0,
              'mseed': {'record_length': 512, 'encoding': 'STEIM2',
                        'filesize': 512, 'dataquality': 'D',
                        'number_of_records': 1, 'byteorder': '>'}}
     stats['starttime'] = UTCDateTime(2000, 1, 1)
     st = Stream([Trace(data=data, header=stats)])
     # Write it.
     st.write(tempfile, format="MSEED")
     # Read it again and delete the temporary file.
     stream = read(tempfile)
     os.remove(tempfile)
     stream.verify()
     # Loop over the attributes to be able to assert them because a
     # dictionary is not a stats dictionary.
     # This also assures that there are no additional keys.
     for key in stats.keys():
         self.assertEqual(stats[key], stream[0].stats[key])

def selectUniqueTraces(tr,args):

    # Test on orizontal component, since if only vertical component exists, 
    # no xcorr on horiz allowed --> crash
  
    st = Stream()
    List = []
    ST = []
    CleanList = []

    for i in range(len(tr)):
        if(tr[i].stats.channel[2:3] == "N"):
           List.append(tr[i].stats.station)

    for i in range(len(tr)):
       a = List.count(tr[i].stats.station)
       if(a > 1):
           ST.append(tr[i].stats.station)

    d = Counter(ST)
    for key in d:
        CleanList.append(key) 
    
    for i in range(len(tr)):
        if CleanList.count(tr[i].stats.station) == 0:
           st.append(tr[i]) 

    return st 

    def write_adjoint_traces(self, path, syn, dat, channel):
        """ Computes adjoint traces from observed and synthetic traces
        """
        nt, dt, _ = self.get_time_scheme(syn)
        nr, _ = self.get_network_size(syn)

        Del = np.loadtxt(path +'/'+ '../../delta_syn_ij')
        rsd = np.loadtxt(path +'/'+ '../../rsd_ij')

        # initialize trace arrays
        adj = Stream()
        for i in range(nr):
            adj.append(Trace(
                data=np.zeros(nt, dtype='float32'),
                header=syn[i].stats))

        # generate adjoint traces
        for i in range(nr):
            for j in range(i):
                si = syn[i].data
                sj = syn[j].data

                adj[i].data += rsd[i,j] * \
                               self.adjoint_dd(si, sj, +Del[i,j], nt, dt)
                adj[j].data -= rsd[i,j] * \
                               self.adjoint_dd(sj, si, -Del[i,j], nt, dt)


        # optional weighting
        adj = self.apply_weights(adj)

        # write adjoint traces
        self.writer(adj, path, channel)

 def test_writeStreamViaObsPy(self):
     """
     Write streams, i.e. multiple files via L{obspy.Trace}
     """
     testdata = np.array([111, 111, 111, 111, 111, 109, 106, 103, 103,
                          110, 121, 132, 139])
     testfile = NamedTemporaryFile().name
     self.file = os.path.join(self.path, '3cssan.reg.8.1.RNON.wav')
     tr = read(self.file, format='WAV')[0]
     np.testing.assert_array_equal(tr.data[:13], testdata)
     # write
     st2 = Stream([Trace(), Trace()])
     st2[0].data = tr.data.copy()       # copy the data
     st2[1].data = tr.data.copy() // 2  # be sure data are different
     st2.write(testfile, format='WAV', framerate=7000)
     # read without giving the WAV format option
     base, ext = os.path.splitext(testfile)
     testfile0 = "%s%03d%s" % (base, 0, ext)
     testfile1 = "%s%03d%s" % (base, 1, ext)
     tr30 = read(testfile0)[0]
     tr31 = read(testfile1)[0]
     self.assertEqual(tr30.stats, tr.stats)
     self.assertEqual(tr31.stats, tr.stats)
     np.testing.assert_array_equal(tr30.data[:13], testdata)
     np.testing.assert_array_equal(tr31.data[:13], testdata // 2)
     os.remove(testfile)
     os.remove(testfile0)
     os.remove(testfile1)

def merge_single(nch,dstart,dend):
  '''Merges traces of one channel to larger traces. Used for cross-correlation'''

  # here you load all the functions you need to use
  from obspy.seg2.seg2 import readSEG2
  from obspy.core import Stream


  # directories:
  dataDir2 = "/import/neptun-radler/STEINACH_feb/"
  dataDir = "/import/three-data/hadzii/STEINACH/STEINACH_longtime/"


  tr = []

  for k in range(dstart, dend, 1):
    fname = '%d' %(k)
    fileName = fname + ".dat" 
    st = readSEG2(dataDir + fileName)
    tr.append(st[nch-1])


  
  new_stream = Stream(traces=tr)
  new_stream.merge(method=1, fill_value='interpolate')


  return new_stream

 def test_writeAndReadDifferentEncodings(self):
     """
     Writes and read a file with different encoding via the obspy.core
     methods.
     """
     npts = 1000
     np.random.seed(815)  # make test reproducable
     data = np.random.randn(npts).astype('float64') * 1e3 + .5
     st = Stream([Trace(data=data)])
     # Loop over some record lengths.
     for encoding, value in ENCODINGS.iteritems():
         seed_dtype = value[2]
         tempfile = NamedTemporaryFile().name
         # Write it once with the encoding key and once with the value.
         st[0].data = data.astype(seed_dtype)
         st.verify()
         st.write(tempfile, format="MSEED", encoding=encoding)
         st2 = read(tempfile)
         del st2[0].stats.mseed
         np.testing.assert_array_equal(st[0].data, st2[0].data)
         del st2
         ms = _MSStruct(tempfile)
         ms.read(-1, 1, 1, 0)
         self.assertEqual(ms.msr.contents.encoding, encoding)
         del ms  # for valgrind
         os.remove(tempfile)