Converting Harmonic decomposition from time to depth

[dip16gphy]

Hi!
I didnt find the command in 'rfpy_hk' script to convert the harmonic decomposition from time to depth domain.
Can u tell how it is done in Audet, P. (2015), Layered crustal anisotropy around the San Andreas Fault near Parkfield, California, J. Geophys. Res. Solid Earth, 120, 3527–3543, doi:10.1002/2014JB011821 ?

[paudetseis]

Hi - the "time-to-depth" conversion in Audet (2015) was done using an older version of the code written in FORTRAN. For development purposes, it was left out of the Python version (RfPy). The necessary functions to do this task appear in the module ccp.py (the function rfpy.ccp.raypath(), which calls some of the other functions in the ccp module), but are not implemented in the scripts. This is something that we may add in the future when most other bugs have been worked out, but is not a trivial task.

If you are comfortable enough with the code base using the API (as opposed to the scripts), you could do this by importing the ccp module and calling the functions on any receiver function Trace object. For instance:

# additional imports
import rfpy.ccp as ccp
from obspy import Stream, Trace

# Insert your own processing to extract RFs, bin them, filter, etc.

# Define velocity model as layers (rather crude model below)
dep = np.array([0., 10., 30., 40., 100.]) # in km
vp = np.array([4.0, 5.0, 5.5, 8.0, 8.1]) # in km/s
vs = np.array([2.3, 2.8, 3.1, 4.5, 4.6]) # in km/s

# Get interpolated velocity profile
idep = np.linspace(dep.min(), dep.max(), nz)
ivp = sp.interpolate.interp1d(dep, vp, kind='linear')(idep)
​ivs = sp.interpolate.interp1d(dep, vs, kind='linear')(idep)

# Number of depth nodes for interpolation
nz = 100

# Initialize streams
rf_stream = Stream()

# Get travel times for Ps, Pps and Pss from velocity model
for tr in rfRstream:

    # Initialize temporary array
    rf_tmp = np.zeros(nz)

   # Travel time calculation
    ttps, ttpps, ttpss, _, _ = ccp.raypath(tr, dep, vp, vs)

   # Get amplitude at each time shift
    for iz, tt in enumerate(ttps):
        amp, _ = ccp.timeshift(tr, tt)
        rf_tmp[iz] = amp

        rf_trace = Trace(data=rf_tmp, header=tr.stats)
        rf_trace.stats.delta = dz
        rf_trace.stats.npts = nz
        rf_trace.stats.taxis = idep

    # Store back into Stream
    rf_stream.append(rf_trace)

# At the end of the loop, you have a new Stream object that 
# contains the time-to-depth converted RF data. You can then 
# use it to do the harmonic decomposition.