def read_enhanced_sac_pz(filename):
    zeros, poles, constant, comments = pz.read_sac_zpk(filename=filename, get_comments=True)
    d = {}
    for line in comments:
        toks = line.split(':', 1)
        if len(toks) == 2:
            temp = toks[0].strip('* \t')
            for k in ('network', 'station', 'location', 'channel', 'start', 'end', 
                      'latitude', 'longitude', 'depth', 'elevation', 'dip', 'azimuth',
                      'input unit', 'output unit'):
                if temp.lower().startswith(k):
                    d[k] = toks[1].strip()

    response = trace.PoleZeroResponse(zeros, poles, constant)

    try:
        channel = Channel(
            nslc=(d['network'], d['station'], d['location'], d['channel']),
            tmin=util.str_to_time(d['start'], format='%Y-%m-%dT%H:%M:%S'),
            tmax=util.str_to_time(d['end'], format='%Y-%m-%dT%H:%M:%S'),
            lat=float(d['latitude']),
            lon=float(d['longitude']),
            elevation=float(d['elevation']),
            depth=float(d['depth']),
            dip=float(d['dip']),
            azimuth=float(d['azimuth']),
            input_unit=d['input unit'],
            output_unit=d['output unit'],
            response=response)
    except:
        raise EnhancedSacPzError('cannot get all required information from file %s' % filename)

    return channel

    def testUSGS(self):

        def is_the_haiti_event(ev):
            assert near(ev.magnitude, 7.0, 0.1)
            assert near(ev.lat, 18.443, 0.01)
            assert near(ev.lon, -72.571, 0.01)
            assert near(ev.depth, 13000., 1.)

        cat = catalog.USGS()

        tmin = util.str_to_time('2010-01-12 21:50:00')
        tmax = util.str_to_time('2010-01-13 03:17:00')

        names = cat.get_event_names(time_range=(tmin, tmax), magmin=5.)
        assert len(names) == 13
        for name in names:
            ev = cat.get_event(name)
            if ev.magnitude >= 7.:
                is_the_haiti_event(ev)
                ident = ev.name

        assert ident is not None
        cat.flush()
        ev = cat.get_event(ident)
        is_the_haiti_event(ev)

def read_data(fn):
    """Read a buletin from the IG CAS website.

    Creates objects of type ResidualMarker (subclass of pyrocko.gui_util.PhaseMarker) and pyrocko.model.Event.
    """
    picks = []
    events = []
    with open(fn, 'r') as f:
        for line in f.readlines():
            odate, eventid, otime, lon, _, lat, _, _, depth, _, _, mag, stat, phase_id, polarity, date, t, a = line.split()
            otime = util.str_to_time('%s %s'%(odate, otime))
            event = model.Event(lat=float(lat), 
                                    lon=float(lon),
                                    depth=float(depth),
                                    time=otime)
            t_pick = util.str_to_time('%s %s'%(date, t))
            pick = ResidualMarker(event=event,
                                  tmin=float(t_pick),
                                  tmax=float(t_pick), 
                                  nslc_ids=[('', stat, '', '*'),],
                                  phasename=phase_id,
                                  polarity=int(polarity),
                                  event_time=otime,
                                  event_hash=event.get_hash())
            
            if abs(t_pick-otime)>30:
                continue
            picks.append(pick)
            events.append(event)

    return picks, events

def readandsplit(infile, box, stations):
   '''
    Load phase files and store data in a myStation.
    Phasename is retrieved from file name: *_*_Phasename, 
    where * can be anything.
    :param infile:      name of the file which is to be read as string.
    :param box:         instance of myStationBox, containing myStation instances.
    :param stations:    File as used with snuffler to import station information.
   '''
   refevent = str_to_time("2010-04-11 22:08:15.500")
   phasename = infile.split('_')[2]
   phasename = phasename.split('.')[0]

   phases = open(infile,'r').readlines()[1::]
   for phase in phases:
       date,time,length,NSLC = phase.split()
       arrival = str_to_time(date+' '+time)-refevent
       netw,stat,loc,comp = NSLC.split('.')
       
       # check if card exists, and if it does, add phase to that card:
       if(box.stationInBox(netw,stat)):
           tmpCard = box.getStationCardForNetworkStation(netw,stat)
           tmpCard.setPhase(phasename,arrival)                               
       # if not, create a new card:
       else:
           newCard=myStation(network=netw,station=stat)
           newCard.setPhase(phasename,arrival)
           newCard.setLatLonEleFromStation(stations)
           box.addStationCard(newCard)

    def testGlobalCMT(self):

        def is_the_haiti_event(ev):
            assert near(ev.magnitude, 7.0, 0.1)
            assert near(ev.lat, 18.61, 0.01)
            assert near(ev.lon, -72.62, 0.01)
            assert near(ev.depth, 12000., 1.)
            assert ev.region.lower() == 'haiti region'

        cat = catalog.GlobalCMT()

        tmin = util.str_to_time('2010-01-12 21:50:00')
        tmax = util.str_to_time('2010-01-13 03:17:00')

        names = cat.get_event_names(time_range=(tmin, tmax), magmin=5.)
        ident = None
        for name in names:
            ev = cat.get_event(name)
            if ev.magnitude > 7:
                is_the_haiti_event(ev)
                ident = ev.name

        assert ident is not None
        cat.flush()
        ev = cat.get_event(ident)
        is_the_haiti_event(ev)

    def testTimeError(self):
        ok = False
        try:
            util.str_to_time('abc')
        except util.TimeStrError:
            ok = True

        assert ok

 def testGeofonMT(self):
     cat = catalog.Geofon()
     tmin = util.str_to_time('2014-01-01 00:00:00')
     tmax = util.str_to_time('2017-01-01 00:00:00')
     events = cat.get_events((tmin, tmax), magmin=8)
     self.assertEqual(len(events), 2)
     mt1, mt2 = [ev.moment_tensor for ev in events]
     angle = moment_tensor.kagan_angle(mt1, mt2)
     self.assertEqual(round(angle - 7.7, 1), 0.0)

def dummy_aware_str_to_time(s, time_format='%Y-%m-%dT%H:%M:%S'):
    try:
        util.str_to_time(s, format=time_format)
    except util.TimeStrError:
        year = int(s[:4])
        if year > this_year + 100:
            return None  # StationXML contained a dummy end date

        raise

 def testTime(self):
     
     for fmt, accu in zip(
         [ '%Y-%m-%d %H:%M:%S.3FRAC', '%Y-%m-%d %H:%M:%S.2FRAC', '%Y-%m-%d %H:%M:%S.1FRAC', '%Y-%m-%d %H:%M:%S' ],
         [ 0.001, 0.01, 0.1, 1.] ):
     
         ta = util.str_to_time('1960-01-01 10:10:10')
         tb = util.str_to_time('2020-01-01 10:10:10')
         
         for i in xrange(10000):
             t1 = ta + random() * (tb-ta)
             s = util.time_to_str(t1, format=fmt)
             t2 = util.str_to_time(s, format=fmt)
             assert abs( t1 - t2 ) < accu

    def testIterTimes(self):

        tmin = util.str_to_time('1999-03-20 20:10:10')
        tmax = util.str_to_time('2001-05-20 10:00:05')

        ii = 0
        for ymin, ymax in util.iter_years(tmin, tmax):
            for mmin, mmax in util.iter_months(ymin, ymax):
                ii += 1
                s1 = util.time_to_str(mmin)
                s2 = util.time_to_str(mmax)

        assert ii == 12*3
        assert s1 == '2001-12-01 00:00:00.000'
        assert s2 == '2002-01-01 00:00:00.000'

    def __init__(self):
        # Set up receiver configuration.

        tab = '''
        HH1  58.500 12.5000  0
        HH2  48.500 12.5000  0
        HH3  48.500  3.5000  0
        HH4  58.500  3.5000  0
        '''.strip()

        receivers = []
        for line_tab in tab.split('\n'):
            station, lat, lon, depth = line_tab.split()
            r = receiver.Receiver(lat, lon, components='neu', name='.%s.' % station)
            receivers.append(r)

        stations = receivers_to_stations(receivers)
        model.dump_stations(stations, 'reference_stations.txt')

        # Composition of the source
        self.olat, self.olon = 52.0000, 9.00000
        self.otime = util.str_to_time('1986-08-22 07:00:00')

        # The gfdb can be chosen within snuffler.
        # This refers to the 'add_parameter' method.
        db = gfdb.Gfdb('fomostos/local1/local1')

        seis = seismosizer.Seismosizer(hosts=['localhost'])
        seis.set_database(db)
        seis.set_effective_dt(db.dt)
        seis.set_local_interpolation('bilinear')
        seis.set_receivers(receivers)
        seis.set_source_location(self.olat, self.olon, self.otime)
        seis.set_source_constraints(0, 0, 0, 0, 0, -1)
        self.seis = seis

    def test_evalresp(self, plot=False):

        resp_fpath = common.test_data_file('test2.resp')

        freqs = num.logspace(num.log10(0.001), num.log10(10.), num=1000)

        transfer = evalresp.evalresp(
            sta_list='BSEG',
            cha_list='BHZ',
            net_code='GR',
            locid='',
            instant=util.str_to_time('2012-01-01 00:00:00'),
            freqs=freqs,
            units='DIS',
            file=resp_fpath,
            rtype='CS')[0][4]

        pz_fpath = common.test_data_file('test2.sacpz')

        zeros, poles, constant = pz.read_sac_zpk(pz_fpath)

        resp = trace.PoleZeroResponse(zeros, poles, constant)

        transfer2 = resp.evaluate(freqs)

        if plot:
            plot_tfs(freqs, [transfer, transfer2])

        assert numeq(transfer, transfer2, 1e-4)

def pdate(s):
    if s.startswith('2599') or s.startswith('2999'):
        return None
    elif s.lower().startswith('no'):
        return None
    else:
        return util.str_to_time(s, format='%Y,%j,%H:%M:%S.OPTFRAC')

    def test_evalresp(self, plot=False):

        testdir = os.path.dirname(__file__)

        freqs = num.logspace(num.log10(0.001), num.log10(10.), num=1000)

        transfer = evalresp.evalresp(sta_list='BSEG',
                          cha_list='BHZ',
                          net_code='GR',
                          locid='',
                          instant=util.str_to_time('2012-01-01 00:00:00'),
                          freqs=freqs,
                          units='DIS',
                          file=os.path.join(testdir, 'response', 'RESP.GR.BSEG..BHZ'),
                          rtype='CS')[0][4]

        pzfn = 'SAC_PZs_GR_BSEG_BHZ__2008.254.00.00.00.0000_2599.365.23.59.59.99999'

        zeros, poles, constant = pz.read_sac_zpk(filename=os.path.join(
            testdir, 'response', pzfn))
        
        resp = trace.PoleZeroResponse(zeros, poles, constant)

        transfer2 = resp.evaluate(freqs)

        if plot:
            import pylab as lab
            lab.plot(freqs, num.imag(transfer))
            lab.plot(freqs, num.imag(transfer2))
            lab.gca().loglog() 
            lab.show()

        assert numeq(transfer, transfer2, 1e-4)

def window(lat, lon, lat_source, lon_source, depth_source, nsta, store): 

        center = Array_center(lat, lon)
        source_reciever_dis = orthodrome.distance_accurate50m_numpy(
            lat_source, lon_source, lat, lon)
    
        t_p = np.array([
            store.t("first(p|P)", (depth_source, int(source_reciever_dis[i])))
            for i in range(0, nsta)])
        t_s = np.array([
           store.t("first(s|S)", (depth_source, int(source_reciever_dis[i])))
         for i in range(0, nsta)])
       
        t_origin = util.str_to_time('2008-02-17 11:06:01.10')
       
        def win_(t_, t_l , t_r, center):
            wind_i = t_origin + t_ - t_l
            wind_e = t_origin + t_ + t_r
            t_o = - t_ + t_[center]
            return wind_i, wind_e, t_o
        
        P_wind_i, P_wind_e, t_op = win_(t_p, 5.0 , 20.0, center)
        S_wind_i, S_wind_e, t_os = win_(t_s, 2.0 , 18.0, center)
        
        return P_wind_i, P_wind_e, t_op , center#, S_wind_i, S_wind_e, t_os

def iload_fh(f, time_format='%Y-%m-%dT%H:%M:%S'):
    zeros, poles, constant, comments = pz.read_sac_zpk(file=f,
                                                       get_comments=True)
    d = {}
    for line in comments:
        toks = line.split(':', 1)
        if len(toks) == 2:
            temp = toks[0].strip('* \t')
            for k in ('network', 'station', 'location', 'channel', 'start',
                      'end', 'latitude', 'longitude', 'depth', 'elevation',
                      'dip', 'azimuth', 'input unit', 'output unit'):

                if temp.lower().startswith(k):
                    d[k] = toks[1].strip()

    response = trace.PoleZeroResponse(zeros, poles, constant)

    try:
        yield EnhancedSacPzResponse(
            codes=(d['network'], d['station'], d['location'], d['channel']),
            tmin=util.str_to_time(d['start'], format=time_format),
            tmax=dummy_aware_str_to_time(d['end']),
            lat=float(d['latitude']),
            lon=float(d['longitude']),
            elevation=float(d['elevation']),
            depth=float(d['depth']),
            dip=float(d['dip']),
            azimuth=float(d['azimuth']),
            input_unit=d['input unit'],
            output_unit=d['output unit'],
            response=response)
    except KeyError as e:
        raise EnhancedSacPzError(
            'cannot get all required information "%s"' % e.args[0])

    def runParallel(inmodel):
        
        print 'does not work properly. EXIT'
        sys.exit(0)
        ProBar = progressbar.ProgressBar(maxval=iterations).start()
        misfits={}
        misfits['pMF']=[]; misfits['sMF']=[]
        misfits['ScsMF']=[]; misfits['ScssMF']=[]

        loadmod=cake.load_model(inmodel)
     
        for latindx, lat in enumerate(_lats):
            for lonindx, lon in enumerate(_lons):
                for zindex, z in enumerate(_depths):
                    #iteration+=1
                    # Start prozess with one event (depth), and one model:
                    eve=model.Event(lat,lon,str_to_time("2010-04-11 22:08:15.500"), 
                       "Spain_Durcal" , z,6.3)
                    [ttpdiff, ttsdiff, ttScsdiff, ttScssdiff] = depthfinder.startup(loadmod, eve, maxdist) 
                   
                    pMF, sMF, ScsMF, ScssMF= map(
                               lambda x: calculateMisfit(x,maxdist), [ttpdiff,ttsdiff,ttScsdiff,ttScssdiff])
                                               
                    resultArray[latindx][lonindx][zindex] = [pMF, sMF, ScsMF, ScssMF]
                    # update progressbar
                    ProBar.update(iteration)
                    identifierstring = inmodel+'.%s.%s.%s'%(lat, lon, z)
                   
                    results[identifierstring]=misfits
        
        try:
            output = open('results.p','w')
            pickle.dump(results, output)
        finally:
            output.close()

    def from_attributes(vals):
        if len(vals) == 14:
            nbasicvals = 7
        else:
            nbasicvals = 4
        nslc_ids, tmin, tmax, kind = Marker.parse_attributes(
            vals[1:1+nbasicvals])

        i = 8
        if len(vals) == 14:
            i = 11

        event_hash = str_to_str_or_none(vals[i-3])
        event_sdate = str_to_str_or_none(vals[i-2])
        event_stime = str_to_str_or_none(vals[i-1])

        if event_sdate is not None and event_stime is not None:
            event_time = util.str_to_time(event_sdate + ' ' + event_stime)
        else:
            event_time = None

        phasename = str_to_str_or_none(vals[i])
        polarity = str_to_int_or_none(vals[i+1])
        automatic = str_to_bool(vals[i+2])
        marker = PhaseMarker(nslc_ids, tmin, tmax, kind, event=None,
                             event_hash=event_hash, event_time=event_time,
                             phasename=phasename, polarity=polarity,
                             automatic=automatic)
        return marker

    def call(self):
         
        self.cleanup()
        
        # Set up receiver configuration.
        tab = '''
        HH  3. 3. 0
        '''.strip()

        receivers = []
        station, lat, lon, depth = tab.split()

        d_north=self.d_north
        d_east=self.d_east
        origin_lat, origin_lon = orthodrome.ne_to_latlon_alternative_method(float(lat), float(lon), d_north, d_east)
        r = receiver.Receiver(lat,lon, components='neu', name='.%s.' % station)
        receivers.append(r)

        # Composition of the source
        otime = util.str_to_time('2000-01-1 00:00:00')
        db = self.db

        seis = seismosizer.Seismosizer(hosts=['localhost'])
        seis.set_database(db)
        seis.set_effective_dt(db.dt)
        seis.set_local_interpolation('bilinear')
        seis.set_receivers(receivers)
        seis.set_source_location( origin_lat, origin_lon, otime)
        seis.set_source_constraints (0, 0, 0, 0 ,0 ,-1)
        self.seis = seis        
        seis = None

        risetime=3; moment=1.
        s = source.Source('bilateral',
        sourceparams_str='0 0 0 %g %g %g %g %g 0 0 0 0 1 %g' % (self.source_depth, moment, self.strike, self.dip, self.rake, risetime))
        self.seis.set_source(s)
        recs = self.seis.get_receivers_snapshot( which_seismograms = ('syn',), which_spectra=(), which_processing='tapered')
        
        trs = []
        for rec in recs:
            if self.save_mseed is True:
                rec.save_traces_mseed(filename_tmpl='%(whichset)s_%(network)s_%(station)s_%(location)s_%(channel)s.mseed' )
            trs.extend(rec.get_traces())
        self.add_traces(trs)
        # Define fade in and out, band pass filter and cut off fader for the TF.
        tfade = self.tfade
        freqlimit = (0.005,.006,1,1.2)
        cut_off_fading = 50
        ntraces = []
        
        for tr in trs:
            TF = STS2()
            
            # Save synthetic trace after transfer function was applied.
            trace_filtered = tr.transfer(tfade, freqlimit, TF, cut_off_fading)            
            # Set new codes to the filtered trace to make it identifiable.
            rename={'e':'BHE','n':'BHN','u':'BHZ'}
            trace_filtered.set_codes(channel=rename[trace_filtered.channel], network='STS2', station='HH', location='syn')
            ntraces.append(trace_filtered)            

    def plot():
        '''
        PLOT PARAMETERS CAN BE ADJUSTED IN THIS METHOD
        '''
        # initialize plotter-instance with pickled stationbox:
        stationBoxFile = open('stationBox.p','r')
        pplotter = plotter(pickle.load(stationBoxFile))
        cmt_lat = 37.10
        cmt_lon = -3.69

        # plot picks (Event coordinates needed to get relative distance):
        CMTevent=model.Event(cmt_lat, cmt_lon, str_to_time("2010-04-11 22:08:15.500"), 
           "Spain_Durcal" , 625500 ,6.3)
        pplotter.plotpicks(CMTevent)

        # plot misfits:
        try:
            pickledResults = open('numpy_results.p','r')
            resultArray = pickle.load(pickledResults)
        finally:
            pickledResults.close()
        
        #pplotter.plotMisfits(data, _depths)
        #.............................................................
        # plot contour misfits for each model:
        # 1st: load misfit array:
        for phaseKey in ['p', 's', 'Scs', 'Scss']:    
            pplotter.plotContourMisfit(resultArray, _lats, _lons, _depths, phaseKey, _models[3])
        

        #.............................................................
        # set time shift file:
        tShiftFiles = glob.glob('map/t_shift*.txt')
        
        if tShiftFiles==[] or (tShiftFiles!=[] and (os.stat(tShiftFiles[0])[8] < os.stat('numpy_results.p')[8])): 
            for phase in ['p','s']:
                for testmodel in _models:
                    print '(re)setting t-shift files for model: {0} and phase:{1}'.format(testmodel, phase)
                    depthfinder.setTshiftFile(CMTevent, np.searchsorted(_lats, CMTevent.lat, ), 
                                np.searchsorted(_lons, CMTevent.lon), resultArray, testmodel, phase, _depths, maxdist )               
        else:
            pass

        for shiftmodel in _models:
            pplotter.plotMisfitShiftMap('s',shiftmodel)

        #.............................................................
        #write horizontal layers to image files;
        #pplotter.saveMisfitLayersAsImages(resultArray, _lats, _lons, _depths, 'models/premlocal1_Stich.nd','p', 630000, 600000 )
        #asdf 

        CMTeventIndexLongitude = getIndexOfValueInArray(_lons, cmt_lon)
        CMTeventIndexLatitude = getIndexOfValueInArray(_lats, cmt_lat)
        pplotter.plotMisfitsForLocation(_models[0], CMTeventIndexLatitude, CMTeventIndexLongitude, resultArray, _depths)


        #.............................................................
        # create map
        subprocess.call("map/durcal_stations.sh&", shell=True)