Python Stream.detrend方法代码示例

# 需要导入模块: from obspy.core.stream import Stream [as 别名]
# 或者: from obspy.core.stream.Stream import detrend [as 别名]
  def section_plot(self, assoc_id, files, seconds_ahead = 5, record_length = 100, channel = 'Z'):
    
    station=self.assoc_db.query(Candidate.sta).filter(Candidate.assoc_id==assoc_id).all()
    sta_list=[]
    for sta, in station:
      sta_list.append(str(sta))
    station_single = self.assoc_db.query(Pick.sta).filter(Pick.assoc_id==assoc_id).filter(Pick.locate_flag == None).all()
    for sta, in station_single:
      sta_list.append(str(sta))
    #print sta_list
      
    eve=self.assoc_db.query(Associated).filter(Associated.id==assoc_id).first()
    # Earthquakes' epicenter
    eq_lat = eve.latitude
    eq_lon = eve.longitude
      
    # Reading the waveforms
    ST = Stream()
    for file in files:
      st = read(file)
      ST += st


    # in case of some seismometer use channel code like BH1, BH2 or BH3, resign the channel code as:
    if channel=='E' or channel=='e':
      Chan='E1'
    elif channel=='N' or channel=='n':
      Chan='N2'
    elif channel=='Z' or channel=='z':
      Chan='Z3'
    else:
      print('Please input component E, e, N, n, Z, or z, the default is Z')
    
    # Calculating distance from headers lat/lon
    ST_new = Stream()#;print ST
    for tr in ST:
      if tr.stats.channel[2] in Chan and tr.stats.station in sta_list:
        if tr.stats.starttime.datetime < eve.ot and tr.stats.endtime.datetime > eve.ot:
          tr.trim(UTCDateTime(eve.ot-timedelta(seconds=seconds_ahead)), UTCDateTime(eve.ot+timedelta(seconds=record_length)))
          ST_new+=tr
    #print ST_new.__str__(extended=True)


    while True:
      ST_new_sta=[]
      for tr in ST_new:
        ST_new_sta.append(tr.stats.station)
      duplicate=list(set([tr for tr in ST_new_sta if ST_new_sta.count(tr)>1]))
      if not duplicate:
        break
      index=[i for (i,j) in enumerate(ST_new_sta) if j==duplicate[-1]]
      i=0
      while True:
        if ST_new[index[i]].stats.npts<ST_new[index[i+1]].stats.npts:
          del ST_new[index[i]]
          break
        elif ST_new[index[i]].stats.npts>=ST_new[index[i+1]].stats.npts:
          del ST_new[index[i+1]]
          break
    #print ST_new.__str__(extended=True)     


    ST_new.detrend('demean')
#     ST_new.filter('bandpass', freqmin=0.1, freqmax=100)

    factor=10
    numRows=len(ST_new)
    segs=[];ticklocs=[];sta=[];circle_x=[];circle_y=[];segs_picks=[];ticklocs_picks=[]
    for tr in ST_new:
      dmax=tr.data.max()
      dmin=tr.data.min()
      data=tr.data/(dmax-dmin)*factor
      t=np.arange(0,round(tr.stats.npts/tr.stats.sampling_rate/tr.stats.delta))*tr.stats.delta # due to the float point arithmetic issue, can not use "t=np.arange(0,tr.stats.npts/tr.stats.sampling_rate,tr.stats.delta)"
      segs.append(np.hstack((data[:,np.newaxis],t[:,np.newaxis])))
      lon,lat = self.tt_stations_db_3D.query(Station3D.longitude,Station3D.latitude).filter(Station3D.sta==tr.stats.station).first()
      distance = int(gps2DistAzimuth(lat,lon,eq_lat,eq_lon)[0]/1000.)  #gps2DistAzimuth return in meters, convert to km by /1000
#       distance=self.assoc_db.query(Candidate.d_km).filter(Candidate.assoc_id==assoc_id).filter(Candidate.sta==tr.stats.station).first()[0]#;print distance,tr.stats.station
      ticklocs.append(distance)
      sta.append(tr.stats.station)
      # DOT plot where picks are picked, notice that for vertical trace plot p is queried from Pick table, s from PickModified table
      # horizontal trace plot p and s queried from PickModified table
      if channel=='Z3':
        picks_p=self.assoc_db.query(Pick.time).filter(Pick.assoc_id==assoc_id).filter(Pick.sta==tr.stats.station).filter(Pick.chan==tr.stats.channel).filter(Pick.phase=='P').all()
        if not picks_p:
          picks_p=self.assoc_db.query(PickModified.time).filter(PickModified.assoc_id==assoc_id).filter(PickModified.sta==tr.stats.station).filter(PickModified.phase=='P').all()
        picks_s=self.assoc_db.query(PickModified.time).filter(PickModified.assoc_id==assoc_id).filter(PickModified.sta==tr.stats.station).filter(PickModified.phase=='S').all()
#         print picks_p,picks_s
      else:
        picks_p=self.assoc_db.query(PickModified.time).filter(PickModified.assoc_id==assoc_id).filter(PickModified.sta==tr.stats.station).filter(PickModified.phase=='P').all()
        picks_s=self.assoc_db.query(PickModified.time).filter(PickModified.assoc_id==assoc_id).filter(PickModified.sta==tr.stats.station).filter(PickModified.phase=='S').all()
#         print picks_p,picks_s
      picks=picks_p+picks_s
#       picks=self.assoc_db.query(PickModified.time).filter(PickModified.assoc_id==assoc_id).filter(PickModified.sta==tr.stats.station).all()
      for pick, in picks:
        index=int((pick-eve.ot+timedelta(seconds=seconds_ahead)).total_seconds()/tr.stats.delta)#;print pick,eve.ot,index,len(data)
        circle_x.append(distance+data[index])
        circle_y.append(t[index])
        # BAR plot where picks are picked  
        t_picks=np.array([t[index],t[index]])
        data_picks=np.array([data.min(),data.max()])
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core.stream import Stream [as 别名]
# 或者: from obspy.core.stream.Stream import detrend [as 别名]
AC = Stream(traces=[BHE[0],BHN[0],BHZ[0]])
ac = AC.copy()


# **Remove the instrument responses of the instruments from the recordings + convert units**
# - convert Ring Laser recordings to nrad/s units using a conversion factor
# - remove the seismometer response using poles and zeros + convert from velocity to acceleration [nm/s^2] in one step
# - trim the traces to make sure start- and endtimes match for both instruments

# In[3]:

RLAS.detrend(type='linear')
RLAS[0].data = RLAS[0].data * 1/6.3191 * 1e-3

AC.detrend(type='linear')
AC.taper(max_percentage=0.05)


paz_sts2 = {'poles': [(-0.0367429 + 0.036754j), (-0.0367429 - 0.036754j)],
            'sensitivity': 0.944019640,
            'zeros': [0j],
            'gain': 1.0}

AC.simulate(paz_remove=paz_sts2, remove_sensitivity=True)

startaim = max([tr.stats.starttime for tr in (AC + RLAS)])
endtaim = min([tr.stats.endtime for tr in (AC + RLAS)])

AC.trim(startaim, endtaim, nearest_sample=True)
RLAS.trim(startaim, endtaim, nearest_sample=True)