def testReadSeisan(self):
        fpath = common.test_data_file('test.seisan_waveform')
        i = 0
        for tr in io.load(fpath, format='seisan'):
            i += 1

        assert i == 39

    def testWriteRead(self):
        now = time.time()
        n = 10
        deltat = 0.1
        
        networks = [ rn(2) for i in range(5) ]
        
        traces1 = [ trace.Trace(rc(networks), rn(4), rn(2), rn(3), tmin=now+i*deltat*n*2, deltat=deltat, ydata=num.arange(n, dtype=num.int32), mtime=now)
            for i in range(3) ]
            
        tempdir = tempfile.mkdtemp()

        for format in ('mseed', 'sac', 'yaff'):
            fns = io.save(traces1, pjoin(tempdir, '%(network)s_%(station)s_%(location)s_%(channel)s'), format=format)

            for fn in fns:
                assert io.detect_format(fn) == format

            traces2 = []
            for fn in fns:
                traces2.extend(io.load(fn, format='detect'))
                
            for tr in traces1:
                assert tr in traces2, 'failed for format %s' % format
                
            for fn in fns:
                os.remove(fn)

        shutil.rmtree(tempdir)

def restitute_evalresp(tr_fn):
    traces = io.load(tr_fn)
    out_traces = []
    for tr in traces:
        try:
            try:
                evalresp = trace.Evalresp(respfile=evalresps['%s.%s'%(tr.station, tr.channel)],
                                          nslc_id=tr.nslc_id,
                                          target='dis')

            except KeyError:
                print 'skip ', '.'.join(tr.nslc_id[1:])
                continue

            if tr.station=='nkc' or tr.station=='zhc':
                t_taper = 30
                f_taper = (0.05, 0.08, 50., 75.),     # frequency domain taper in [hz]
            else:
                t_taper = 5.
                f_taper = (0.3, 0.6, 50., 75.),     # frequency domain taper in [hz]

            displacement =  tr.transfer(
                t_taper,                       # rise and fall of time domain taper in [s]
                *f_taper,     # frequency domain taper in [hz]
                transfer_function=evalresp, 
                invert=True)

        except trace.TraceTooShort:
            continue
        out_traces.append(displacement)
    tr_fn = tr_fn.replace(inputdir, outputdir)
    print tr_fn
    io.save(out_traces, tr_fn)

    del traces

    def load_data(self, force=False):
        file_changed = False
        if not self.data_loaded or force:
            logger.debug('loading data from file: %s' % self.abspath)
            
            for itr, tr in enumerate(io.load(self.abspath, format=self.format, getdata=True, substitutions=self.substitutions)):
                if itr < len(self.traces):
                    xtr = self.traces[itr]
                    if xtr.mtime != tr.mtime or xtr.tmin != tr.tmin or xtr.tmax != tr.tmax:
                        logger.debug('file may have changed since last access (trace number %i has changed): %s' % (itr, self.abspath))
                        self.remove(xtr)
                        self.traces.remove(xtr)
                        xtr.file = None
                        self.traces.append(tr)
                        self.add(tr)
                        tr.file = self
                        file_changed = True
                    else:
                        xtr.ydata = tr.ydata

                else:
                    self.traces.add(tr)
                    self.add(tr)
                    logger.debug('file may have changed since last access (new trace found): %s' % self.abspath)
                    file_changed = True
            self.data_loaded = True
        return file_changed

    def testReadKan(self):
        fpath = common.test_data_file('01.kan')
        i = 0
        for tr in io.load(fpath, format='kan'):
            i += 1

        assert i == 1

    def call(self):
        """Main work routine of the snuffling."""

        self.cleanup()

        view = self.get_viewer()
        pile = self.get_pile()

        tmin, tmax = view.get_time_range()
        if self.useevent:
            markers = view.selected_markers()
            if len(markers) != 1:
                self.fail("Exactly one marker must be selected.")
            marker = markers[0]
            if not isinstance(marker, EventMarker):
                self.fail("An event marker must be selected.")

            ev = marker.get_event()

            lat, lon = ev.lat, ev.lon
        else:
            lat, lon = self.lat, self.lon

        data = iris_ws.ws_station(
            lat=lat, lon=lon, minradius=self.minradius, maxradius=self.maxradius, timewindow=(tmin, tmax), level="chan"
        )

        stations = iris_ws.grok_station_xml(data, tmin, tmax)
        networks = set([s.network for s in stations])

        t2s = util.time_to_str
        dir = self.tempdir()
        fns = []
        for net in networks:
            nstations = [s for s in stations if s.network == net]
            selection = sorted(iris_ws.data_selection(nstations, tmin, tmax))
            if selection:
                for x in selection:
                    logger.info("Adding data selection: %s.%s.%s.%s %s - %s" % (tuple(x[:4]) + (t2s(x[4]), t2s(x[5]))))

                try:
                    d = iris_ws.ws_bulkdataselect(selection)
                    fn = pjoin(dir, "data-%s.mseed" % net)
                    f = open(fn, "w")
                    f.write(d)
                    f.close()
                    fns.append(fn)

                except urllib2.HTTPError:
                    pass

        all_traces = []
        for fn in fns:
            try:
                traces = list(io.load(fn))

                all_traces.extend(traces)

            except io.FileLoadError, e:
                logger.warning("File load error, %s" % e)

    def testReadGSE1(self):
        fpath = common.test_data_file('test1.gse1')
        i = 0
        for tr in io.load(fpath, format='detect'):
            i += 1

        assert i == 19

    def call(self):
        '''Main work routine of the snuffling.'''
        
        self.cleanup()

        view = self.get_viewer()
        pile = self.get_pile()

        tmin, tmax = view.get_time_range()
        if self.useevent:
            markers = view.selected_markers()
            if len(markers) != 1:
                self.fail('Exactly one marker must be selected.')
            marker = markers[0]
            if not isinstance(marker, EventMarker):
                self.fail('An event marker must be selected.')

            ev = marker.get_event()
            
            lat, lon = ev.lat, ev.lon
        else:
            lat, lon = self.lat, self.lon
        
        print lat, lon, self.minradius, self.maxradius, util.time_to_str(tmin), util.time_to_str(tmax)

        data = iris_ws.ws_station(lat=lat, lon=lon, minradius=self.minradius, maxradius=self.maxradius, 
                                                     timewindow=(tmin,tmax), level='chan' )
        stations = iris_ws.grok_station_xml(data, tmin, tmax)
        networks = set( [ s.network for s in stations ] )
        
        dir = self.tempdir()
        fns = []
        for net in networks:
            nstations = [ s for s in stations if s.network == net ]
            selection = sorted(iris_ws.data_selection( nstations, tmin, tmax ))
            if selection:
                for x in selection:
                    print x
                
                try:
                    d = iris_ws.ws_bulkdataselect(selection)
                    fn = pjoin(dir,'data-%s.mseed' % net) 
                    f = open(fn, 'w')
                    f.write(d)
                    f.close()
                    fns.append(fn)
                except urllib2.HTTPError:
                    pass

        newstations = []
        for sta in stations:
            if not view.has_station(sta):
                print sta
                newstations.append(sta)

        view.add_stations(newstations)
        for fn in fns:
            traces = list(io.load(fn))
            self.add_traces(traces)

    def testReadSEGY(self):
        fpath = common.test_data_file('test2.segy')
        i = 0
        for tr in io.load(fpath, format='segy'):
            assert tr.meta['orfield_num'] == 1111
            i += 1

        assert i == 24

    def testReadCSS(self):
        wfpath = common.test_data_file('test_css1.w')  # noqa
        fpath = common.test_data_file('test_css.wfdisc')
        i = 0
        for tr in io.load(fpath, format='css'):
            i += 1

        assert i == 1

    def testReadGcf(self):
        fpath = common.test_data_file('test.gcf')

        i = 0
        for tr in io.load(fpath, format='gcf'):
            i += 1

        assert i == 1

    def test_ahfull_kiwi(self):
        setup = load(filename=common.test_data_file(
            'test_ahfull_kiwi_setup.yaml'))
        trs_ref = io.load(common.test_data_file(
            'test_ahfull_kiwi_traces.mseed'))

        for i, s in enumerate(setup.setups):
            d3d = math.sqrt(s.x[0]**2 + s.x[1]**2 + s.x[2]**2)

            tlen = d3d / s.vs * 2

            n = int(num.round(tlen / s.deltat))

            out_x = num.zeros(n)
            out_y = num.zeros(n)
            out_z = num.zeros(n)

            ahfullgreen.add_seismogram(
                s.vp, s.vs, s.density, 1000000.0, 1000000.0, s.x, s.f, s.m6,
                'displacement',
                s.deltat, 0.,
                out_x, out_y, out_z,
                ahfullgreen.Gauss(s.tau))

            trs = []
            for out, comp in zip([out_x, out_y, out_z], 'NED'):
                tr = trace.Trace(
                    '', 'S%03i' % i, 'P', comp,
                    deltat=s.deltat, tmin=0.0, ydata=out)

                trs.append(tr)

            trs2 = []

            for cha in 'NED':

                t1 = g(trs, 'S%03i' % i, cha)
                t2 = g(trs_ref, 'S%03i' % i, cha)

                tmin = max(t1.tmin, t2.tmin)
                tmax = min(t1.tmax, t2.tmax)

                t1 = t1.chop(tmin, tmax, inplace=False)
                t2 = t2.chop(tmin, tmax, inplace=False)

                trs2.append(t2)

                d = 2.0 * num.sum((t1.ydata - t2.ydata)**2) / \
                    (num.sum(t1.ydata**2) + num.sum(t2.ydata**2))

                if d >= 0.02:
                    print(d)
                    # trace.snuffle([t1, t2])

                assert d < 0.02

    def testReadSUDS(self):
        fpath = common.test_data_file('test.suds')
        i = 0
        for tr in io.load(fpath, format='detect'):
            i += 1

        assert i == 251

        stations = suds.load_stations(fpath)

        assert len(stations) == 91

    def test_to_obspy_trace(self):
        traces = io.load(common.test_data_file('test1.mseed'))
        for tr in traces:
            obs_tr = tr.to_obspy_trace()

            assert isinstance(obs_tr, obspy.Trace)
            assert obs_tr.data.size == tr.data_len()

            obs_stats = obs_tr.stats
            for attr in ('network', 'station', 'location', 'channel'):
                assert obs_stats.__getattr__(attr) == tr.__getattribute__(attr)

    def load_headers(self, mtime=None):
        logger.debug('loading headers from file: %s' % self.abspath)
        if mtime is None:
            self.mtime = os.stat(self.abspath)[8]
        
        self.remove(self.traces)
        for tr in io.load(self.abspath, format=self.format, getdata=False, substitutions=self.substitutions):
            self.traces.append(tr)
            tr.file = self

        self.add(self.traces)

        self.data_loaded = False
        self.data_use_count = 0

def restitute_pz(tr_fn):
    traces = io.load(tr_fn)
    out_traces = []
    for tr in traces:
        try:
            try:
                zeros, poles, constant = pz.read_sac_zpk(pole_zeros['%s.%s'%(tr.station, tr.channel)])

            except keyerror:
                print 'skip ', '.'.join(tr.nslc_id[1:])
                continue

            zeros.append(0.0j)

            digitizer_gain = 1e6    
            constant *= digitizer_gain

            # for the conversion of hz-> iw:
            nzeros = len(zeros)
            npoles = len(poles)
            constant *= (2*num.pi)**(npoles-nzeros)

            if tr.station=='nkc' or tr.station=='zhc':
                constant *= normalization_factors[tr.station]
                t_taper = 30
                f_taper = (0.05, 0.08, 50., 75.),     # frequency domain taper in [hz]
            else:
                t_taper = 5.
                f_taper = (0.3, 0.6, 50., 75.),     # frequency domain taper in [hz]

            print tr.station, constant
            pz_transfer = trace.polezeroresponse(zeros, poles, constant)
            displacement =  tr.transfer(
                t_taper,                       # rise and fall of time domain taper in [s]
                *f_taper,     # frequency domain taper in [hz]
                transfer_function=pz_transfer, 
                invert=true)

        except trace.tracetooshort:
            continue
        out_traces.append(displacement)
    tr_fn = tr_fn.replace(inputdir, outputdir)
    print tr_fn
    io.save(out_traces, tr_fn)

    del traces

 def get_traces_pyrocko(self, x,z):
     fns = self.dump_traces(x,z, format='mseed')
     traces = []
     for igm, fn in enumerate(fns):
         ig = igm+1
         if fn:
             for tr in io.load(fn):
                 ix = 1 + int(round((x-self.firstx)/self.dx))
                 iz = 1 + int(round((z-self.firstz)/self.dz))
                 gridx = self.firstx + (ix-1)*self.dx
                 gridz = self.firstz + (iz-1)*self.dz
                 sx = util.base36encode(ix)
                 sz = util.base36encode(iz)
                 tr.meta = {'x':gridx, 'z':gridz, 'ig':ig}
                 tr.set_codes(network=sz, station=sx, channel='%i' % ig)
                 traces.append(tr)
                 
     return traces

# write the incoming data stream to 'traces.mseed'
with open('traces.mseed', 'wb') as file:
    file.write(request_waveform.read())

# request meta data
request_response = fdsn.station(
    site='geofon', selection=selection, level='response')

# save the response in YAML and StationXML format
request_response.dump(filename='responses_geofon.yaml')
request_response.dump_xml(filename='responses_geofon.xml')

# Loop through retrieved waveforms and request meta information
# for each trace
traces = io.load('traces.mseed')
displacement = []
for tr in traces:
    polezero_response = request_response.get_pyrocko_response(
        nslc=tr.nslc_id,
        timespan=(tr.tmin, tr.tmax),
        fake_input_units='M')
    # *fake_input_units*: required for consistent responses throughout entire
    # data set

    # deconvolve transfer function
    restituted = tr.transfer(
        tfade=2.,
        freqlimits=(0.01, 0.1, 1., 2.),
        transfer_function=polezero_response,
        invert=True)

from pyrocko import io

traces = io.load('test.mseed')
t = traces[0]
print 'original:', t

# extract a copy of a part of t
extracted = t.chop(t.tmin+10, t.tmax-10, inplace=False)
print 'extracted:', extracted

# in-place operation modifies t itself
t.chop(t.tmin+10, t.tmax-10)
print 'modified:', t

    def download(self, event, directory='array_data', timing=None, length=None,
                 want='all', force=False, prefix=False, dump_config=False,
                 get_responses=False):
        """:param want: either 'all' or ID as string or list of IDs as strings
        """
        use = []
        #ts = {}
        unit = 'M'
        if all([timing, length]) is None:
            raise Exception('Define one of "timing" and "length"')
        prefix = prefix or ''
        directory = pjoin(prefix, directory)
        if not os.path.isdir(directory):
            os.mkdir(directory)
        pzresponses = {}
        logger.info('download data: %s at %sN %sE' % (
            event.name, event.lat, event.lon))
        for site, array_data_provder in self.providers.items():
            logger.info('requesting data from site %s' % site)
            for array_id, codes in array_data_provder.items():
                if array_id not in want and want != ['all']:
                    continue
                sub_directory = pjoin(directory, array_id)
                logger.info("%s" % array_id)
                codes = array_data_provder[array_id]
                if not isinstance(codes, list):
                    codes = [codes]
                selection = [
                    c + tuple((event.time, event.time+1000.)) for c in codes]
                logger.debug('selection: %s' % selection)
                try:
                #    if site=='bgr':
                #        st = ws.station(url='http://eida.bgr.de/', selection=selection)
                #    else:
                #        st = ws.station(site=site, selection=selection)
                    st = ws.station(site=site, selection=selection)
                except ws.EmptyResult as e:
                    logging.error('No results: %s %s. skip' % (e, array_id))
                    continue
                except ValueError as e:
                    logger.error(e)
                    logger.error('...skipping...')
                    continue

                stations = st.get_pyrocko_stations()
                min_dist = min(
                    [ortho.distance_accurate50m(s, event) for s in stations])
                max_dist = max(
                    [ortho.distance_accurate50m(s, event) for s in stations])

                mod = cake.load_model(crust2_profile=(event.lat, event.lon))
                if length:
                    tstart = 0.
                    tend = length
                elif timing:
                    tstart = timing[0].t(mod, (event.depth, min_dist))
                    tend = timing[1].t(mod, (event.depth, max_dist))
                selection = [
                    c + tuple((event.time + tstart, event.time + tend)
                              ) for c in codes]
                try:
                    d = ws.dataselect(site=site, selection=selection)
                    store.remake_dir(sub_directory, force)
                    store.remake_dir(pjoin(sub_directory, 'responses'), force)
                    fn = pjoin(sub_directory, 'traces.mseed')
                    with open(fn, 'w') as f:
                        f.write(d.read())
                        f.close()
                    if get_responses:
                        trs = io.load(fn, getdata=False)
                        logger.info('Request responses from %s' % site)
                        if progressbar:
                            pb = progressbar.ProgressBar(maxval=len(trs)).start()
                        for i_tr, tr in enumerate(trs):
                            try:
                                st = ws.station(
                                    site=site, selection=selection, level='response')
                                pzresponse = st.get_pyrocko_response(
                                    nslc=tr.nslc_id,
                                    timespan=(tr.tmin, tr.tmax),
                                    fake_input_units=unit)
                                pzresponse.regularize()
                            except fdsnstation.NoResponseInformation as e:
                                logger.warn("no response information: %s" % e)
                                pzresponse = None
                                pass
                            except fdsnstation.MultipleResponseInformation as e:
                                logger.warn("MultipleResponseInformation: %s" % e)
                                pzresponse = None
                                pass
                            pzresponses[tr.nslc_id] = pzresponse
                            pzresponses[tr.nslc_id].dump(filename=pjoin(
                                sub_directory,
                                'responses',
                                'resp_%s.yaml' % '.'.join(tr.nslc_id)))
                            if progressbar:
                                pb.update(i_tr)
                        if progressbar:
                            pb.finish()
                    model.dump_stations(
#.........这里部分代码省略.........