Preparing RETREAT body waves for Receiver-Function Estimation

Work performed on an IntelMac laptop -- computer codes and scripts will be posted in successive blog posts. This post documents the preparation of body wave data from the RETREAT experiment. You can use the concepts to guide the preparation of data for other stations.

The directory Teleseismic has subdirectories PIIR, VLC, VOLR, etc, one for each station name. The station subdirectories have data files in SAC format, but without the source location in the header. There are three steps in preparing a station's data for processing. First, we associate events with data-file triplets. The extraction by RDSEED from the SEED files should have synchronized the data files for each component of motion. This first step can be tedious, however, because there will be more events than datafile triplets whenever the station was down, and more data files than events when the data files are fragmented. The goal is to prune the files "events" to only the events that have data, and the file "infiles1" to be the triplets of datafile names -- infiles1 will have three times as many lines and events, natch.

cd PIIR (from the Teleseismic directory)
ls *.SAC > infiles_raw
cp infiles_raw infiles1
cp ../events.julian events
te events &
te infiles1 &

events_julian is a master file of events used to request data from the IRIS DMC for the experiment, based on a JWeed event file:

2003 356 MHDF 2003 12 22 19:15:56.000 35.706 -121.102 7.6 MB 6.1 MS 6.4
2003 359 MHDF 2003 12 25 07:11:11.590 8.416 -82.824 33 MB 6.0 MS 6.4
2003 359 MHDF 2003 12 25 20:42:33.720 -22.252 169.488 10 MB 6.3 MS 6.3
2003 360 MHDF 2003 12 26 01:56:52.440 28.995 58.311 10 MB 6.0 MS 6.8
2003 360 MHDF 2003 12 26 21:26:04.100 -22.273 169.314 10 MB 6.1 MS 6.8
2003 361 MHDF 2003 12 27 16:00:59.450 -22.015 169.766 10 MB 6.1 MS 7.1
2003 361 MHDF 2003 12 27 22:38:01.880 -21.672 169.835 10 MB 5.8 MS 6.7

Note that the month/day/year format of the JWeed event file is replaced by year julianday, which matches the output filename format of rdseed.


The user now opens infiles1 and events in a text editor and synchronizes their lines. It is smart to keep an Xterminal open to run SAC or pql. Sometimes there will be more than three datafiles associated with an event start time, and you must display their contents in order to choose the three files to use, or else to merge some files.

pql yyyy.ddd.* &

is a good command for this. Or

echo "read yyyy.ddd.*; qdp 2000; p1" | sac

because you can retrieve previous commands in a Unix shell and replace the yyyy.ddd for different cases.

The second step is to run the code "rotate_sac_puteventinheader" or "rotate_sac." this code reads the info in infiles1 and events. It generates a SAC macro file that will read the data triplets, put the earthquake hypocentre into the SAC headers (NOT for sac_rotate, only for sac_rotate_puteventinheader), cut the files to 600-s length and rotate the horizontals to radial and transverse direction. You run

rotate_sac_puteventinheader
/bin/mv outfile sac_rotate
echo exit >> sac_rotate
sac < sac_rotate

and the files *yyyy.ddd.hh.*BH[ENZ]*.SAC are converted to truncated and rotated data files "yyyy.ddd.hh.mm.bh[zrt]" Now comes the third step of processing. The code examine_sac reads the data files and composes a SAC macro to visualize the data in raw and bandpassed traces, in order to assess the signal to noise ratio of the P wave.

ls *bh? >! infiles2
echo 3 >! c_c_c
examine_sac < c_c_c
/bin/mv outfile1 in_recfunk
/bin/mv outfile sac_examine
sed -e 's/plot1;pause/plotpk markall on; writehdr over/' sac_examine >! sac_examine_pick
sac < sac_examine_pick

The SAC-macro sac_examine displays data from each event in raw, lowpassed and highpassed traces. Oftentimes a Pwave is evident in some bandpasses and not others -- I typically keep events where the low or high bandpass has SNR>2, but SNR.le.1 elsewhere. The last "sed" command converts this display macro into a phase-picking macro. "sac_examine_pick" asks you to click "A" on the desired start of a RF data window (the duration of the window is user-input in the RF codes, but the start time of analysis is most convenient if stored in the SAC header). The SAC macro marks "A" in the lowpassed traces at the point 120s from the start of the record, the typical timing of the P wave. A PKP or PP wave will arrive later than this, except for events that are so far away that JWEED neglected to use the Pdiff as a reference. GCARC~150 is the transition, but Im not certain where exactly this transition occurs. For P and Pdiff selection, you can leave the cursor close to 120s after trace onset, and adjust to capture the P wave with a 5-10 seconds to spare. For the PKP phase, you must move the cursor to later in the record.

The data-evaluation and marking exercise is done with the A and Q keys.
If the P or PKP wave is determined "analyzable" click A to mark the start time of the data window, then Q to move to the next record. SAC marks all data traces with the A-time. If you make a mistake in selecting A, either too early or too late, just mark it again.
If the P or PKP wave is weak or bad (e.g. the PP overlaps the PKP too much), then simply click Q and SAC will leave the datafile unmarked. This flags the event and datafiles to be discarded. If you click Q by accident before clicking A before a good P phase, you must return to the data later by running the relevant lines of sac_examine_pick. Write down the event time (yyyy.ddd.hh) and go back to it later.
examine_sac takes the data files and orders them in increasing epicentral distance. This way the user will assess the data first from nearby earthquakes with high-frequency P, moving on to low-frequency Pdiff, to an interval where only the largest quakes have visible Pdiff, and then to GCARC~120 where the PKP phase precedes the PP phase by enough time to be useful for RF analysis. Finally the time-reference for JWeed shifts and the data interval has PKP/PKIKP starting at 120s.

After the desired P waves have been identified with A-markers in their SAC headers, it is time to prune the list of events down to the chosen few. To prune the in_recfunk file, 

grep bh in_recfunk >! in_recpick_raw
echo stop >> in_recpick_raw
recfunk_pick_through 

--> reads in_recpick_raw and writes in_recpick, a list of wildcarded filenames that includes only those whose headers define a data window to analyze. Then this in_recpick can be used by recfunk_pick to compute simple RF sweeps (no migration, only epicentral and back-azimuths bin-sums).

recfunk_pick

There are GMT scripts in the Teleseismic directory, some files from recfunk_pick can be plotted to see what the station RFs are like:

../Plotss outr_baz.grid outt_baz.grid 0.75 plot_baz
cp plot.ps plot_baz.ps

../Plot_epi3 outr_epi.grid outt_epi.grid 0.75 plot_epi
cp plot.ps plot_epi.ps