Python obspy.Trace方法代碼示例

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def flex_cut_trace(trace, cut_starttime, cut_endtime, dynamic_npts=0):
    """
    not cut strictly(but also based on the original trace length)

    :param trace: input trace
    :type trace: obspy.Trace
    :param cut_starttime: starttime of cutting
    :type cut_starttime: obspy.UTCDateTime
    :param cut_endtime: endtime of cutting
    :type cut_endtime: obspy.UTCDateTime
    """
    if not isinstance(trace, Trace):
        raise TypeError("flex_cut_trace method only accepts obspy.Trace "
                        "as the first argument")

    delta = trace.stats.delta
    cut_starttime = cut_starttime - dynamic_npts * delta
    cut_endtime = cut_endtime + dynamic_npts * delta
    trace.trim(cut_starttime, cut_endtime) 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def convert_adj_to_trace(adj):
    """
    Convert AdjointSource to Trace,for internal use only
    """
    meta = {}

    tr = Trace()
    tr.data = deepcopy(adj.adjoint_source)
    tr.stats.starttime = adj.starttime
    tr.stats.delta = adj.dt

    tr.stats.channel = adj.component
    tr.stats.station = adj.station
    tr.stats.network = adj.network
    tr.stats.location = adj.location

    meta["adj_src_type"] = adj.adj_src_type
    meta["misfit"] = adj.misfit
    meta["min_period"] = adj.min_period
    meta["max_period"] = adj.max_period

    return tr, meta 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def convert_trace_to_adj(tr, meta):
    """
    Convert Trace to AdjointSource, for internal use only, with
    meta data information
    """
    minp = meta["min_period"]
    maxp = meta["max_period"]
    adj_src_type = meta["adj_src_type"]
    misfit = meta["misfit"]

    dt = tr.stats.delta
    component = tr.stats.channel
    adj = AdjointSource(adj_src_type, misfit, dt, minp, maxp, component)

    adj.adjoint_source = deepcopy(tr.data)
    adj.station = tr.stats.station
    adj.network = tr.stats.network
    adj.location = tr.stats.location
    adj.starttime = tr.stats.starttime

    return adj 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def tw2utc(tw, trace):
    """Convert time window to UTC time window

    :param tw: tuple of two values, both can be a string (see :func:`time2utc`)
        or a list of strings in which case the latest starttime and earliest
        endtime is taken.
    :param trace: Trace object with stats entries
    """
    starttime = None
    for val in tw:
        if isinstance(val, (list, tuple)):
            times = [time2utc(v, trace, starttime=starttime) for v in val]
            t = max(times) if starttime is None else min(times)
        else:
            t = time2utc(val, trace, starttime=starttime)
        if starttime is None:
            starttime = t
    return starttime, t 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def check_sample(stream):
    """
    Returns sampling rate of traces in stream.

    :type stream:`~obspy.core.stream.Stream` object. 
    :param stream: Stream containing one or more day-long trace 
    :return: List of sampling rates in stream

    """
    if type(stream) == obspy.core.trace.Trace:
        return stream
    else:
        freqs = []	
        for tr in stream:
            freqs.append(tr.stats.sampling_rate)

    freq = max(set(freqs),key=freqs.count)
    for tr in stream:
        if tr.stats.sampling_rate != freq:
            stream.remove(tr)

    return stream 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def downsample(stream,freq):
    """ 
    Downsamples stream to specified samplerate.

    Uses Obspy.core.trace.decimate if mod(sampling_rate) == 0. 
    :type stream:`~obspy.core.trace.Stream` or `~obspy.core.trace.Trace` object.
    :type freq: float
    :param freq: Frequency to which waveforms in stream are downsampled
    :return: Downsampled trace or stream object
    :rtype: `~obspy.core.trace.Trace` or `~obspy.core.trace.Trace` object.
    """
    
    # get sampling rate 
    if type(stream) == obspy.core.stream.Stream:
        sampling_rate = stream[0].stats.sampling_rate
    elif type(stream) == obspy.core.trace.Trace:
        sampling_rate = stream.stats.sampling_rate

    if sampling_rate == freq:
        pass
    else:
        stream.interpolate(freq,method="weighted_average_slopes")	

    return stream 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def _upsample(self, trace, upfactor):
        """
        Upsample a data stream by a given factor, prior to decimation. The
        upsampling is done using a linear interpolation.

        Parameters
        ----------
        trace : obspy Trace object
            Trace to be upsampled
        upfactor : int
            Factor by which to upsample the data in trace

        Returns
        -------
        out : obpsy Trace object
            Upsampled trace

        """

        data = trace.data
        dnew = np.zeros(len(data) * upfactor - (upfactor - 1))
        dnew[::upfactor] = data
        for i in range(1, upfactor):
            dnew[i::upfactor] = float(i) / upfactor * data[1:] \
                         + float(upfactor - i) / upfactor * data[:-1]

        out = Trace()
        out.data = dnew
        out.stats = trace.stats
        out.stats.npts = len(out.data)
        out.stats.starttime = trace.stats.starttime
        out.stats.sampling_rate = int(upfactor * trace.stats.sampling_rate)

        return out 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def filter_trace(tr, pre_filt):
    """
    Perform a frequency domain taper mimicing the behavior during the
    response removal, without a actual response removal.

    :param tr: input trace
    :param pre_filt: frequency array(Hz) in ascending order, to define
        the four corners of filter, for example, [0.01, 0.1, 0.2, 0.5].
    :type pre_filt: Numpy.array or list
    :return: filtered trace
    """
    if not isinstance(tr, Trace):
        raise TypeError("First Argument should be trace: %s" % type(tr))
    if len(pre_filt) != 4:
        raise ValueError("Length of filter must be 4(corner frequencies)")
    if not check_array_order(pre_filt, order="ascending"):
        raise ValueError("Frequency band should be in ascending order: %s"
                         % pre_filt)

    data = tr.data.astype(np.float64)
    origin_len = len(data)
    if origin_len == 0:
        return

    # smart calculation of nfft dodging large primes
    nfft = _npts2nfft(len(data))

    fy = 1.0 / (tr.stats.delta * 2.0)
    freqs = np.linspace(0, fy, nfft // 2 + 1)

    # Transform data to Frequency domain
    data = np.fft.rfft(data, n=nfft)
    data *= c_sac_taper(freqs, flimit=pre_filt)
    data[-1] = abs(data[-1]) + 0.0j
    # transform data back into the time domain
    data = np.fft.irfft(data)[0:origin_len]
    # assign processed data and store processing information
    tr.data = data 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def plot_adjoint_source(adjsrc, win_times=None,
                        obs_tr=None, syn_tr=None,
                        figname=None):
    """
    Plot adjoint source for multiple windows

    :param figname: output figure file name
    :type figname: str
    :param adjsrc: adjoint source
    :type adjsrc: pyadjoint.AdjointSource
    :param adjsrc
    :return:
    """
    if not isinstance(adjsrc, AdjointSource):
        raise ValueError("Input adjsrc should be type of "
                         "pyadjoint.AdjointSource")

    if obs_tr is None or syn_tr is None:
        plot_only_adjoint(adjsrc, win_times)
    else:
        if not isinstance(obs_tr, Trace):
            raise ValueError("Input obs_tr should be type of obspy.Trace")
        if not isinstance(syn_tr, Trace):
            raise ValueError("Input syn_tr should be type of obspy.Trace")
        if win_times is None:
            raise ValueError("Input win_tims should be specified as time "
                             "of windows")
        plot_adjoint_and_data(adjsrc, win_times, obs_tr, syn_tr)

    if figname is None:
        plt.show()
    else:
        plt.savefig(figname) 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def observed_energy(stream, rho, df, coda_normalization=None, fs=4, tolerance=1):
    """
    Return trace with total spectral energy density of three component stream

    :param stream: stream of a 3 component seismogram
    :param rho: density (kg/m**3)
    :param df: filter width in Hz
    :param fs: free surface correction (default: 4)
    :param tolerance: the number of samples the length of the traces
        in the 3 component stream may differ (default: 1)
    :return: trace with total energy density"""
    data = [energy1c(tr.data, rho, df, fs=fs) for tr in stream]
    Ns = [len(d) for d in data]
    if max(Ns) - min(Ns) > tolerance:
        msg = ('traces for one stream have different lengths %s. Tolerance '
               ' is %d samples') % (Ns, tolerance)
        raise SkipError(msg)
    elif max(Ns) - min(Ns) > 0:
        data = [d[:min(Ns)] for d in data]
    data = np.sum(data, axis=0)
    tr = obspy.Trace(data=data, header=stream[0].stats)
    tr.stats.channel = tr.stats.channel[:2] + 'X'
    if coda_normalization is not None:
        sl = tr.slice(tr.stats.origintime + coda_normalization[0],
                      tr.stats.origintime + coda_normalization[1])
        tr.data = tr.data / np.mean(sl.data)
    return tr 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def match_trace(trace,stream):
    """
    Matches trace in stream that begin at the same time UTC. 

    Removes matched trace from stream for faster matching.
    
    :type trace:`~obspy.core.trace.Trace` object. 
    :param trace: Day-long trace 
    :type stream:`~obspy.core.stream.Stream` object. 
    :param stream: Stream containing one or more day-long trace 
    :Returns: trace from stream object that has same starting time 
    :rtype:`~obspy.core.trace.Trace` object. 
     """
    
    # max time difference between starting sample 0 minutes
    max_time = trace.stats.delta 

    matched_trace = False

    for ii,tr in enumerate(stream):
        if np.abs(tr.stats.starttime - trace.stats.starttime) <= max_time and \
        len(tr.data) == len(trace.data):
            matched_trace = tr
            stream.pop(ii)
            break

    return matched_trace,stream 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def remove_resp(arr,stats,inv):
    """
    Removes instrument response of cross-correlation

    :type arr: numpy.ndarray 
    :type stats: `~obspy.core.trace.Stats` object.
    :type inv: `~obspy.core.inventory.inventory.Inventory`
    :param inv: StationXML file containing response information
    :returns: cross-correlation with response removed
    """	
    
    def arr_to_trace(arr,stats):
        tr = obspy.Trace(arr)
        tr.stats = stats
        tr.stats.npts = len(tr.data)
        return tr

    # prefilter and remove response
    
    st = obspy.Stream()
    if len(arr.shape) == 2:
        for row in arr:
            tr = arr_to_trace(row,stats)
            st += tr
    else:
        tr = arr_to_trace(arr,stats)
        st += tr
    min_freq = 1/tr.stats.npts*tr.stats.sampling_rate
    min_freq = np.max([min_freq,0.005])
    pre_filt = [min_freq,min_freq*1.5, 0.9*tr.stats.sampling_rate, 0.95*tr.stats.sampling_rate]
    st.attach_response(inv)
    st.remove_response(output="VEL",pre_filt=pre_filt) 

    if len(st) > 1: 
        data = []
        for tr in st:
            data.append(tr.data)
        data = np.array(data)
    else:
        data = st[0].data
    return data 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def preprocess(trace,percent=0.01,max_len=20.):   
    """
    Removes linear trend and mean, normalizes and tapers Obspy trace. 
    
    :type trace:`~obspy.core.trace.Trace` object.   
    :type: percent: float, optional
    :param percent: percent window on each end of trace to taper
    :return: Processed trace 
    """
    trace.detrend(type='constant')
    trace.detrend(type='simple')
    percent = trace.stats.sampling_rate * 20 / trace.stats.npts
    trace.taper(max_percentage=percent,max_length=max_len) 	

    return trace 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def _downsample(self, stream, sr, upfactor=None):
        """
        Downsample the stream to the specified sampling rate.

        Parameters
        ----------
        stream : obspy Stream object
            Contains list of Trace objects to be downsampled

        sr : int
            Output sampling rate

        Returns
        -------
        stream : obspy Stream object
            Contains list of Trace objects, with Traces downsampled / resampled
            where necessary and possible.

        """

        for trace in stream:
            if sr != trace.stats.sampling_rate:
                if (trace.stats.sampling_rate % sr) == 0:
                    trace.filter("lowpass", freq=float(sr) / 2.000001,
                                 corners=2, zerophase=True)
                    trace.decimate(factor=int(trace.stats.sampling_rate / sr),
                                   strict_length=False,
                                   no_filter=True)
                elif self.resample and upfactor is not None:
                    # Check the upsampled sampling rate can be decimated to sr
                    if int(trace.stats.sampling_rate * upfactor) % sr != 0:
                        raise util.BadUpfactorException
                    stream.remove(trace)
                    trace = self._upsample(trace, upfactor)
                    trace.filter("lowpass", freq=float(sr) / 2.000001,
                                 corners=2, zerophase=True)
                    trace.decimate(factor=int(trace.stats.sampling_rate / sr),
                                   strict_length=False,
                                   no_filter=True)
                    stream += trace
                else:
                    msg = "Mismatched sampling rates - cannot decimate data - "
                    msg += "to resample data, set .resample = True and choose"
                    msg += " a suitable upfactor"
                    print(msg)

        return stream 

# 需要導入模塊: import obspy [as 別名]
# 或者: from obspy import Trace [as 別名]
def calculate_misfit(_tr1, _tr2, taper_flag=True, taper_percentage=0.05,
                     correlation_flag=True):
    """
    Calculate the misfit between two traces
    :param tr1: trace 1
    :type tr1: Obspy.Trace
    :param tr2: trace 2
    :type tr2: Obspy.Trace
    :param taper_flag: taper the seismogram or not
    :type taper_flag: bool
    :param taper_percentage: the taper percentage
    :type taper_percentage: float
    """
    if not isinstance(_tr1, Trace):
        raise TypeError("Input tr1(type:%s) must be type of obspy.Trace"
                        % type(_tr1))
    if not isinstance(_tr2, Trace):
        raise TypeError("Input tr2(type:%s) must be type of obspy.Trace"
                        % type(_tr2))

    tr1 = _tr1.copy()
    tr2 = _tr2.copy()

    starttime = max(tr1.stats.starttime, tr2.stats.starttime)
    endtime = min(tr1.stats.endtime, tr2.stats.endtime)
    sampling_rate = min(tr1.stats.sampling_rate, tr2.stats.sampling_rate)
    npts = int((endtime - starttime) * sampling_rate)

    tr1.interpolate(sampling_rate, starttime=starttime, npts=npts)
    tr2.interpolate(sampling_rate, starttime=starttime, npts=npts)

    if taper_flag:
        tr1.taper(max_percentage=taper_percentage, type='hann')
        tr2.taper(max_percentage=taper_percentage, type='hann')

    corr_min = cross_correlation(tr1.data, tr2.data)
    err_max = least_squre_error(tr1.data, tr2.data)

    # coverage
    tr1_cover = trace_length(tr1) / trace_length(_tr1)
    tr2_cover = trace_length(tr2) / trace_length(_tr2)

    # amplitude diff
    twdiff = [i / sampling_rate for i in range(npts)]
    amp_ref = np.sum(np.abs(tr1.data) + np.abs(tr2.data)) / (2 * npts)
    wdiff = (tr1.data - tr2.data) / amp_ref

    return {"tr1_coverage": tr1_cover, "tr2_coverage": tr2_cover,
            "correlation": corr_min, "error": err_max,
            "time_array": twdiff, "diff_array": wdiff}