Python Stream.sort方法代码示例

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
def axisem2mseed_all(path):
    """
    change .dat files into MSEED format
    """
    global test_param

    t = UTCDateTime(0)
    traces = []

    for file in glob.iglob(os.path.join(path, '*.dat')):
        stationID = file.split('/')[-1].split('_')[0]
        chan = file.split('/')[-1].split('_')[-1].split('.')[0]
        dat = np.loadtxt(file)
        npts = len(dat[:,0])
        stats = {'network': 'SG',
                 'station': stationID,
                 'location': '',
                 'channel': chan,
                 'npts': npts,
                 'sampling_rate': (npts - 1.)/(dat[-1,0] - dat[0,0]),
                 'starttime': t + dat[0,0],
                 'mseed' : {'dataquality': 'D'}}
        traces.append(Trace(data=dat[:,1], header=stats))

    st = Stream(traces)
    st.sort()
    fname =  os.path.join(path, 'seismograms.mseed')
    print fname
    st.write(fname, format='MSEED')

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
num_stations = 0
exceptions = []
for station in STATIONS:
    try:
        # we request 60s more at start and end and cut them off later to avoid
        # a false trigger due to the tapering during instrument correction
        tmp = client.waveform.getWaveform(NET, station, "", CHANNEL, T1 - 180,
                                          T2 + 180, getPAZ=True,
                                          getCoordinates=True)
    except Exception, e:
        exceptions.append("%s: %s" % (e.__class__.__name__, e))
        continue
    st.extend(tmp)
    num_stations += 1
st.merge(-1)
st.sort()

summary = []
summary.append("#" * 79)
summary.append("######## %s  ---  %s ########" % (T1, T2))
summary.append("#" * 79)
summary.append(st.__str__(extended=True))
if exceptions:
    summary.append("#" * 33 + " Exceptions  " + "#" * 33)
    summary += exceptions
summary.append("#" * 79)

trig = []
mutt = []
if st:
    # preprocessing, backup original data for plotting at end

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
def preprocess(db, stations, comps, goal_day, params, responses=None):

    datafiles = {}
    output = Stream()
    for station in stations:
        datafiles[station] = {}
        net, sta = station.split('.')
        gd = datetime.datetime.strptime(goal_day, '%Y-%m-%d')
        files = get_data_availability(
            db, net=net, sta=sta, starttime=gd, endtime=gd)
        for comp in comps:
            datafiles[station][comp] = []
        for file in files:
            if file.comp[-1] not in comps:
                continue
            fullpath = os.path.join(file.path, file.file)
            datafiles[station][file.comp[-1]].append(fullpath)

    for istation, station in enumerate(stations):
        net, sta = station.split(".")
        for comp in comps:
            files = eval("datafiles['%s']['%s']" % (station, comp))
            if len(files) != 0:
                logging.debug("%s.%s Reading %i Files" %
                              (station, comp, len(files)))
                stream = Stream()
                for file in sorted(files):
                    st = read(file, dytpe=np.float,
                              starttime=UTCDateTime(gd),
                              endtime=UTCDateTime(gd)+86400)
                    tmp = st.select(network=net, station=sta, component=comp)
                    if not len(tmp):
                        for tr in st:
                            tr.stats.network = net
                        st = st.select(network=net, station=sta, component=comp)
                    else:
                        st = tmp
                    for tr in st:
                        tr.data = tr.data.astype(np.float)
                    stream += st
                    del st
                stream.sort()
                stream.merge(method=1, interpolation_samples=3, fill_value=None)
                stream = stream.split()

                logging.debug("Checking sample alignment")
                for i, trace in enumerate(stream):
                    stream[i] = check_and_phase_shift(trace)

                logging.debug("Checking Gaps")
                if len(getGaps(stream)) > 0:
                    max_gap = 10
                    only_too_long = False
                    while getGaps(stream) and not only_too_long:
                        too_long = 0
                        gaps = getGaps(stream)
                        for gap in gaps:
                            if int(gap[-1]) <= max_gap:
                                stream[gap[0]] = stream[gap[0]].__add__(stream[gap[1]], method=1,
                                                                        fill_value="interpolate")
                                stream.remove(stream[gap[1]])
                                break
                            else:
                                too_long += 1
                        if too_long == len(gaps):
                            only_too_long = True
                stream = stream.split()
                taper_length = 20.0  # seconds
                for trace in stream:
                    if trace.stats.npts < 4 * taper_length * trace.stats.sampling_rate:
                        stream.remove(trace)
                    else:
                        trace.detrend(type="demean")
                        trace.detrend(type="linear")
                        trace.taper(max_percentage=None, max_length=1.0)

                if not len(stream):
                    logging.debug(" has only too small traces, skipping...")
                    continue

                for trace in stream:
                    logging.debug(
                        "%s.%s Highpass at %.2f Hz" % (station, comp, params.preprocess_highpass))
                    trace.filter("highpass", freq=params.preprocess_highpass, zerophase=True)

                    if trace.stats.sampling_rate != params.goal_sampling_rate:
                        logging.debug(
                            "%s.%s Lowpass at %.2f Hz" % (station, comp, params.preprocess_lowpass))
                        trace.filter("lowpass", freq=params.preprocess_lowpass, zerophase=True, corners=8)

                        if params.resampling_method == "Resample":
                            logging.debug("%s.%s Downsample to %.1f Hz" %
                                          (station, comp, params.goal_sampling_rate))
                            trace.data = resample(
                                trace.data, params.goal_sampling_rate / trace.stats.sampling_rate, 'sinc_fastest')

                        elif params.resampling_method == "Decimate":
                            decimation_factor = trace.stats.sampling_rate / params.goal_sampling_rate
                            if not int(decimation_factor) == decimation_factor:
                                logging.warning("%s.%s CANNOT be decimated by an integer factor, consider using Resample or Lanczos methods"
#.........这里部分代码省略.........

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
    # write as ASCII file (encoding=0)
    #outfile_name += ".mseed"
    #st.write(outfile_name, format='MSEED', encoding=0, reclen=256)
    st.write(outfile_name, format='MSEED', encoding=11, reclen=256, byteorder='>')

    # Show that it worked, convert NumPy character array back to string
    st1 = read(outfile_name)
    print(st1[0].data.tostring())

    exit(0)

# 2) Construct Stream object
st = read(infile_paths[0])
for i in range(1, len(infile_paths)):
    st += read(infile_paths[i])
st = st.sort()

st[0].stats['mseed'] = {'dataquality': 'D'}
st[0].data = np.int32(st[0].data)

# 3) Convert
if type_name == "SLIST" or type_name == "slist":
    outfile_name += ".slist"
    st.write(outfile_name, format="SLIST")
    exit(0)
if type_name == "TSPAIR" or type_name == "tspair":
    outfile_name += ".tspair"
    st.write(outfile_name, format="TSPAIR")
    exit(0)
if type_name == "MSEED" or type_name == "mseed":
    outfile_name += ".mseed"

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
import sm_utils


#Need to add comments to this code

debug = True

files = glob.glob('data/*')

st = Stream()
for curfile in files:
	st+=read(curfile)

#Make everything overlap and nice
st._cleanup()
st.sort(['network','station','channel','starttime'])


pers = np.linspace(0.01, 20, 300)
for idx, tr in enumerate(st):
	tr.detrend('demean')
	tr.taper(max_percentage=0.05, type='cosine')
	tr.data = tr.data
	tr.filter('highpass',freq=1.0)
	rs = NigamJennings(tr.data,1/tr.stats.sampling_rate,pers,damping=0.5, units="m/s/s")
	spec, ts, acc, vel, dis = rs.evaluate()
	fig = plt.figure(1, figsize=(8,8))
	plt.loglog(spec["Period"],spec["Pseudo-Acceleration"]*0.1/9.81,'k',linewidth=2)
	plt.xlim([.1, 10])
	plt.xlabel('Period (s)', fontsize=18)
	plt.ylabel('Acceleration (g)', fontsize=18)

# 需要导入模块: from obspy.core import Stream [as 别名]
# 或者: from obspy.core.Stream import sort [as 别名]
def preprocess(db, stations, comps, goal_day, params, responses=None):
    """
    Fetches data for each ``stations`` and each ``comps`` using the
    data_availability table in the database.

    To correct for instrument responses, make sure to set ``remove_response``
    to "Y" in the config and to provide the ``responses`` DataFrame.

    :Example:

    >>> from msnoise.api import connect, get_params, preload_instrument_responses
    >>> from msnoise.preprocessing import preprocess
    >>> db = connect()
    >>> params = get_params(db)
    >>> responses = preload_instrument_responses(db)
    >>> st = preprocess(db, ["YA.UV06","YA.UV10"], ["Z",], "2010-09-01", params, responses)
    >>> st
     2 Trace(s) in Stream:
    YA.UV06.00.HHZ | 2010-09-01T00:00:00.000000Z - 2010-09-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples
    YA.UV10.00.HHZ | 2010-09-01T00:00:00.000000Z - 2010-09-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples

    :type db: :class:`sqlalchemy.orm.session.Session`
    :param db: A :class:`~sqlalchemy.orm.session.Session` object, as
        obtained by :func:`msnoise.api.connect`.
    :type stations: list of str
    :param stations: a list of station names, in the format NET.STA.
    :type comps: list of str
    :param comps: a list of component names, in Z,N,E,1,2.
    :type goal_day: str
    :param goal_day: the day of data to load, ISO 8601 format: e.g. 2016-12-31.
    :type params: class
    :param params: an object containing the config parameters, as obtained by
        :func:`msnoise.api.get_params`.
    :type responses: :class:`pandas.DataFrame`
    :param responses: a DataFrame containing the instrument responses, as
        obtained by :func:`msnoise.api.preload_instrument_responses`.
    :rtype: :class:`obspy.core.stream.Stream`
    :return: A Stream object containing all traces.
    """
    datafiles = {}
    output = Stream()
    for station in stations:
        datafiles[station] = {}
        net, sta = station.split('.')
        gd = datetime.datetime.strptime(goal_day, '%Y-%m-%d')
        files = get_data_availability(
            db, net=net, sta=sta, starttime=gd, endtime=gd)
        for comp in comps:
            datafiles[station][comp] = []
        for file in files:
            if file.comp[-1] not in comps:
                continue
            fullpath = os.path.join(file.path, file.file)
            datafiles[station][file.comp[-1]].append(fullpath)

    for istation, station in enumerate(stations):
        net, sta = station.split(".")
        for comp in comps:
            files = eval("datafiles['%s']['%s']" % (station, comp))
            if len(files) != 0:
                logger.debug("%s.%s Reading %i Files" %
                              (station, comp, len(files)))
                stream = Stream()
                for file in sorted(files):
                    try:
                        st = read(file, dytpe=np.float,
                              starttime=UTCDateTime(gd),
                              endtime=UTCDateTime(gd)+86400)
                    except:
                        logger.debug("ERROR reading file %s" % file)
                        continue
                    for tr in st:
                        if len(tr.stats.channel) == 2:
                            tr.stats.channel += tr.stats.location
                            tr.stats.location = "00"
                    tmp = st.select(network=net, station=sta, component=comp)
                    if not len(tmp):
                        for tr in st:
                            tr.stats.network = net
                        st = st.select(network=net, station=sta, component=comp)
                    else:
                        st = tmp
                    for tr in st:
                        tr.data = tr.data.astype(np.float)
                        tr.stats.network = tr.stats.network.upper()
                        tr.stats.station = tr.stats.station.upper()
                        tr.stats.channel = tr.stats.channel.upper()

                    stream += st
                    del st
                stream.sort()
                try:
                    # HACK not super clean... should find a way to prevent the
                    # same trace id with different sps to occur
                    stream.merge(method=1, interpolation_samples=3, fill_value=None)
                except:
                    logger.debug("Error while merging...")
                    traceback.print_exc()
                    continue
                stream = stream.split()
#.........这里部分代码省略.........