Python ASDFDataSet.process方法代码示例

# 需要导入模块: from pyasdf import ASDFDataSet [as 别名]
# 或者: from pyasdf.ASDFDataSet import process [as 别名]
def test_processing_multiprocessing(example_data_set):
    """
    Tests the processing using multiprocessing.
    """
    def null_processing(st, inv):
        return st

    data_set = ASDFDataSet(example_data_set.filename)
    output_filename = os.path.join(example_data_set.tmpdir, "output.h5")
    # Do not actually do anything. Apply an empty function.
    data_set.process(null_processing, output_filename,
                     {"raw_recording": "raw_recording"})

    del data_set
    data_set = ASDFDataSet(example_data_set.filename)
    out_data_set = ASDFDataSet(output_filename)

    assert data_set == out_data_set

# 需要导入模块: from pyasdf import ASDFDataSet [as 别名]
# 或者: from pyasdf.ASDFDataSet import process [as 别名]
        st.taper(max_percentage=0.05, type="hann")

        st.interpolate(sampling_rate=sampling_rate, starttime=starttime,
                       npts=npts)

        station_latitude = inv[0][0].latitude
        station_longitude = inv[0][0].longitude
        _, baz, _ = gps2DistAzimuth(station_latitude, station_longitude,
                                    event_latitude, event_longitude)

        components = [tr.stats.channel[-1] for tr in st]
        if "N" in components and "E" in components:
            st.rotate(method="NE->RT", back_azimuth=baz)

        # Convert to single precision to save space.
        for tr in st:
            tr.data = np.require(tr.data, dtype="float32")

        return st

    tag_name = "preprocessed_%is_to_%is" % (int(min_period), int(max_period))

    tag_map = {
        "raw_recording": tag_name
    }

    ds.process(process_function, tag_name + ".h5", tag_map=tag_map)

# Important when running with MPI as it might otherwise not be able to finish.
del ds

# 需要导入模块: from pyasdf import ASDFDataSet [as 别名]
# 或者: from pyasdf.ASDFDataSet import process [as 别名]
def process_synt(asdf_fn, outputfn, filter_band, old_tag=None, new_tag=None):
    # read in dataset
    ds = ASDFDataSet(asdf_fn)

    max_period = filter_band[1]
    min_period = filter_band[0]
    if min_period > max_period:
        raise ValueError("filter_band incorrect: min_period > max_period")
    f2 = 1.0 / max_period
    f3 = 1.0 / min_period
    f1 = 0.8 * f2
    f4 = 1.2 * f3
    pre_filt = (f1, f2, f3, f4)

    # read in event
    event = ds.events[0]
    origin = event.preferred_origin() or event.origins[0]
    event_latitude = origin.latitude
    event_longitude = origin.longitude
    event_time = origin.time

    # Figure out these parameters somehonw!
    starttime = event_time
    npts = 3600
    sampling_rate = 1.0

    def process_function(st, inv):
        st.detrend("linear")
        st.detrend("demean")
        st.taper(max_percentage=0.05, type="hann")

        # Perform a frequency domain taper like during the response removal
        # just without an actual response...
        for tr in st:
            data = tr.data.astype(np.float64)

            # smart calculation of nfft dodging large primes
            from obspy.signal.util import _npts2nfft
            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:len(data)]
            # assign processed data and store processing information
            tr.data = data

        st.detrend("linear")
        st.detrend("demean")
        st.taper(max_percentage=0.05, type="hann")

        st.interpolate(sampling_rate=sampling_rate, starttime=starttime,
                    npts=npts)

        components = [tr.stats.channel[-1] for tr in st]
        if "N" in components and "E" in components:
            station_latitude = float(inv[0][0].latitude)
            station_longitude = float(inv[0][0].longitude)
            _, baz, _ = gps2DistAzimuth(station_latitude, station_longitude,
                                    event_latitude, event_longitude)

            st.rotate(method="NE->RT", back_azimuth=baz)

        # Convert to single precision to save space.
        for tr in st:
            tr.data = np.require(tr.data, dtype="float32")

        return st

    tag_map = {
            old_tag : new_tag
    }
    
    # process
    ds.process(process_function, outputfn, tag_map=tag_map)

# 需要导入模块: from pyasdf import ASDFDataSet [as 别名]
# 或者: from pyasdf.ASDFDataSet import process [as 别名]
        station_latitude = float(inv[0][0].latitude)
        station_longitude = float(inv[0][0].longitude)
        _, baz, _ = gps2DistAzimuth(station_latitude, station_longitude,
                                    event_latitude, event_longitude)

        components = [tr.stats.channel[-1] for tr in st]
        if "N" in components and "E" in components:
            st.rotate(method="NE->RT", back_azimuth=baz)

        # Convert to single precision to save space.
        for tr in st:
            tr.data = np.require(tr.data, dtype="float32")

        return st

    new_tag = "proc_obsd_%i_%i" % (int(min_period), int(max_period))

    tag_map = {
        old_tag : new_tag 
    }

    outputfn = eventname + "." + new_tag + ".h5"
    outputfn = os.path.join(outputdir, outputfn)
    if os.path.exists(outputfn):
        os.remove(outputfn)
    ds.process(process_function, outputfn, tag_map=tag_map)

t2=time.time()

print "Elapsed time:", t2-t1