text: #!/home/zross/bin/python 
  
  import numpy as np
  import sys
  import glob
  import obspy
  import pylab as plt
  import json
  import random
  
  class Arrival():
      def __init__(self, net=None, sta=None, time=None, phase=None,
          dist=None, resid=None):
          self.net = net
          self.sta = sta
          self.time = time
          self.phase = phase
          self.dist = dist
          self.resid = resid
  
  class Event():
      def __init__(self, arrivals = None):
          if arrivals is not None:
              self.arrivals = arrivals
          else:
              self.arrivals = []
  
  def get_unassociated_trigs(origin_time, triggers, trig_meta):
      t_start = origin_time - obspy.UTCDateTime(0) - 60.0
      t_stop = t_start + 120.
      idx = np.where((triggers >= t_start) & (triggers < t_stop))[0]
      trigs = {}
      for x in idx:
          if trig_meta[x][1] not in trigs:
              trigs[trig_meta[x][1]] = []
          trigs[trig_meta[x][1]].append((trig_meta[x][3], trig_meta[x][4]))
      return trigs
  
  def plot_seismicity(catalog, params):
      import pandas as pd
  
      print('Reading fault file in GMT format, please wait...')
  
      # list to store fault segments
      faults = []
  
      # preallocate to track fault pts within segment
      maxpts = 1600000 # based on number of lines in file
      flats = np.zeros(maxpts)
      flons = np.zeros(maxpts)
      fsegs = np.zeros(maxpts,dtype='int')
      nn = -1
      nseg=-1
  
      # loop over lines
      with open(params['fault_file']) as f:
          for line in f:
  
              # header line that gives number of points in segment
              if line.startswith('Pline'):
                  nseg+=1
  
              # fault point line
              elif line.startswith('-1'):
                  nn+=1
                  lineS = line.split()
                  flons[nn]=float(lineS[0])
                  flats[nn]=float(lineS[1])
                  fsegs[nn]=nseg
  
      # covert to dataframe
      fault_df = pd.DataFrame()
      fault_df['flon']=flons[:nn+1]
      fault_df['flat']=flats[:nn+1]
      fault_df['fseg']=fsegs[:nn+1]
      print('Done, {:} faults read'.format(nseg+1))
  
      from mpl_toolkits.basemap import Basemap, shiftgrid, cm
      fig = plt.figure(figsize=(10,10))
      ax = plt.gca()
      lat0, lat1 = params['lat_min'], params['lat_max']
      clat = (lat0+lat1)/2.
      lon0, lon1 = params['lon_min'], params['lon_max']
      clon = (lon0+lon1)/2.
  
      proj = 'merc'
      epsg = 4269
      m = Basemap(llcrnrlon=lon0,llcrnrlat=lat0,urcrnrlon=lon1,urcrnrlat=lat1,
                  resolution='h',projection=proj,lat_0=clat,lon_0=clon, ax=ax,
                  epsg=epsg)
      m.drawcoastlines()
      m.fillcontinents(color='white', lake_color='paleturquoise')
      m.drawparallels(np.arange(32, 38, 1.), labels=[1,0,0,1])
      m.drawmeridians(np.arange(-120, -114, 1.), labels=[1,0,0,1])
      m.drawmapboundary(fill_color='paleturquoise')
  
  
      xpixels = 5000
      service = 'World_Shaded_Relief'
      #m.arcgisimage(service=service, xpixels = xpixels, verbose= False)
  
      # plot faults
      ifaults = (fault_df.flat >= lat0)&(fault_df.flat <= lat1) & (
                  fault_df.flon >= lon0)&(fault_df.flon <= lon1)
      for g, v in fault_df[ifaults].groupby('fseg'):
          m.plot(v.flon.values,v.flat.values,'-k',lw=1.0,latlon=True)
  
      lon = []
      lat = []
      for event in cat:
          lon.append(event.origins[0].longitude)
          lat.append(event.origins[0].latitude)
      #with open("datasets/cahuilla_sum.nll", 'r') as f:
      #    for line in f:
      #        temp = line.split()
      #        lon.append(float(temp[11]))
      #        lat.append(float(temp[9]))
  
      m.scatter(lon, lat, 0.5, marker='o', color='r', latlon=True, zorder=10)
      stla = []
      stlo = []
      with open(params["station_file"], 'r') as f:
          for line in f:
              temp = line.split()
              stla.append(float(temp[2]))
              stlo.append(float(temp[3]))
      m.scatter(stlo, stla, 50, marker='^', color='blue', latlon=True, zorder=10)
      plt.tight_layout()
      plt.savefig("detection_map.png", dpi=320)
      plt.show()
  
  
  if __name__ == "__main__":
      if len(sys.argv) != 2:
          print("Usage: python phaselink_plot.py config_json")
          print("E.g. python phaselink_plot.py params.json")
          sys.exit()
      with open(sys.argv[1], "r") as f:
          params = json.load(f)
  
      triggers = []
      trig_meta = []
      if params['plot_unassociated']:
          print("Reading unassociated triggers...")
          with open(params['gpd_file'], 'r') as f:
              for line in f:
                  net, sta, phase, time, prob, dur = line.split()
                  if float(prob) < params['pr_min'] or \
                     float(dur) < params['trig_dur_min']:
                      continue
                      trig_type = 0
                  else:
                      trig_type = 1
                  triggers.append(obspy.UTCDateTime(time) - obspy.UTCDateTime(0))
                  trig_meta.append((net, sta, phase, obspy.UTCDateTime(time),
                      trig_type))
          idx = np.argsort(triggers)
          triggers = np.array([triggers[x] for x in idx])
          trig_meta = [trig_meta[x] for x in idx]
  
      print("Now building catalog")
  
      #nll_summary_file = "%s/%s" % \
      #    (params['nlloc_loc_path'], params['nlloc_sum_file'])
      #cat = obspy.io.nlloc.core.read_nlloc_hyp(nll_summary_file)
      nll_files = glob.glob("%s/*.*.*.*.*.hyp" % params['nlloc_loc_path'])
      cat = obspy.Catalog()
      for fname in nll_files:
          try:
              cat += obspy.read_events(fname)
          except:
              continue
      random.shuffle(nll_files)
  
      for event in cat:
          print(event.preferred_origin().time)
      print(cat)
      print()
  
      if params['plot_seismicity']:
          plot_seismicity(cat, params)
  
  
      for fname in nll_files:
          cat = obspy.read_events(fname)
          event = cat[0]
          origin = event.preferred_origin()
          origin_time = origin.time
          print(event)
          print(origin)
  
          if params['plot_unassociated']:
              trigs = get_unassociated_trigs(origin_time, triggers, trig_meta)
  
          # Build id_map for join between arrivals and picks
          picks = {}
          sta_order = []
          dist_count = 0
          for arrival in origin.arrivals:
              pick = arrival.pick_id.get_referred_object()
              sta = pick.waveform_id.station_code
              phase = arrival.phase
              time = pick.time
              #if arrival.distance <= params['dist_cutoff_radius']:
              #    dist_count += 1
              if abs(arrival.time_residual) > params['max_t_resid']:
                  flag = 1
              else:
                  flag = 0
              if sta not in picks:
                  picks[sta] = [(time, phase, flag)]
                  sta_order.append(sta)
              else:
                  picks[sta].append((time, phase, flag))
  
          #if dist_count < params['dist_cutoff_n_min']:
          #    print("Skipping event, only %d phases within radius %.2f" % \
          #        (dist_count, params['dist_cutoff_radius']))
          #    continue
  
          # Plot results
          fig, ax = plt.subplots(1,1,figsize=(30,30))
          colors = {0: 'lime', 1: 'yellow'}
          count = 0
          for sta in sta_order:
  
              st = obspy.read("%s/%04d/%03d/*.%s.*" % \
                  (params['wf_path'], origin_time.year, origin_time.julday, sta),
                  starttime=origin_time-60, endtime=origin_time+60)
              st.detrend()
              st.filter(type='bandpass', freqmin=3.0, freqmax=20)
              for tr in st:
                  ax.plot(np.arange(tr.data.size)*tr.stats.delta, \
                          tr.data/np.max(tr.data) + count, c='k', lw=1)
                  ax.text(125, count, sta)
                  if params['plot_unassociated']:
                      if sta in trigs:
                          for pick, t_type in trigs[sta]:
                              #tr_slice = tr.slice(starttime=pick,
                              #                    endtime=pick+1.0)
                              #ax.plot(np.arange(tr_slice.data.size) \
                              #    * tr.stats.delta + (pick - origin_time) + 60.,
                              #    tr_slice.data/np.max(tr.data) + count,
                              #    c=colors[t_type], lw=1)
                              ax.plot(pick-tr.stats.starttime, 0,
                                      marker="|", c=colors[t_type])
  
                  for pick, phase, flag in picks[sta]:
                      if phase == 'P':
                          color = 'r'
                      else:
                          color = 'b'
                      #if flag:
                      #    color = 'limegreen'
                      ax.plot([pick-tr.stats.starttime, pick-tr.stats.starttime], [count-0.75, count+0.75], c=color)
                  count += 1
          plt.show()
          print()
  