Preprocessing Status

compare.ipynb

process_data.py

import numpy as np
import obspy
from obspy_sdf import SDFDataSet

ds = SDFDataSet("observed_data.h5")

# I don't yet understand how the times/npts are derived.
starttime = obspy.UTCDateTime("2009-01-15T07:27:30.019481Z")
npts = 5708


# Loop over both period sets. This will result in two files. It could also be
# saved to the same file.
for min_period, max_period in ((27.0, 60.0), (60.0, 120.0)):
    print 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)

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

        st.attach_response(inv)
        st.remove_response(output="DISP", pre_filt=pre_filt, zero_mean=False,
                           taper=False)

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

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

        # 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)