Introduction to Seismic Analysis Code (SAC)
by Zhigang Peng
This website contains a brief tutorial on Seismic Analysis Code (SAC).
It is part of the mini course series offered to the incoming Geophysics graduate students at GT.
This tutorial consists of the following sections:
- The current version of SAC is already installed at "/usr/local/geophysics/sac/sac". To run it, type "sac" in the command line.
If you get something like
% sac
SEISMIC ANALYSIS CODE [03/01/2005 (Version 100.00)]
Copyright 1995 Regents of the University of California
SAC>
It it working.
- If it is not working for some reason, check the following two variables.
% echo $SACAUX
/usr/local/geophysics/sac/sac/aux
% echo $PATH
- You should have SACAUX defined as "/usr/local/geophysics/sac/sac/aux", and "/usr/local/geophysics/bin/" in
your path. If not, please add the following two lins in your ${HOME}/.bashrc files to make it work
properly.
# added for Seismic Analysis Code (SAC)
export SACAUX=/usr/local/geophysics/sac/sac/aux
export PATH=/usr/local/geophysics/bin:/usr/local/geophysics/sac/local_bin:/usr/local/geophysics/sac/sac/bin:$PATH
- For users in other institutions that are using tcsh, please use the following examples to add corresponding
lines in your ${HOME}/.tcshrc.
# added for Seismic Analysis Code (SAC)
setenv SACHOME /usr/local/geophysics/sac/sac
setenv SACAUX $SACHOME/aux
set path=($SACHOME/bin $path)
- The standard data format for SAC is binary. A binary SAC file contains a fixed-length ASCII data (header), which describe the
subsequent data in binary format with variable length.
- For the seismic data this means a single data component recorded at a single seismic station.
SAC does not currently work on multiplexed data.
- There is an issue with the SAC binary data. In UNIX, SAC data is stored according to the
"big-endian" byte order (high-order byte of the number is stored in memory at the lowest address).
In LINUX (Intel processors), it is stored according to the "little-endian" byte order (see also Endianness).
- Since most of us are working in LINUX now, it is best to convert the data to the "little-endian" byte order
when you obtain a SAC data that is in UNIX's "big-endian" byte order.
- Use a program called "sacsun2linux" to do the job. Make sure that you only run the program
once.
- Let's start with an example. For GT user, please untar the following file to your home directory,
or any subdirectory where you would like to put the tutorial material.
% tar zxvf /usr/local/geophysics/sac/tutorial/Sac_Tutorial_2009.tar.gz
% cd SAC_Tutorial/SAC_Tutorial_WF
- For outside user, please download the example files from the following link:
% wget http://geophysics.eas.gatech.edu/people/zpeng/Teaching/Sac_Tutorial_2009.tar.gz
% tar zxvf Sac_Tutorial_2009.tar.gz
% cd SAC_Tutorial/SAC_Tutorial_WF
- Next, we will read some seismograms recorded during a temporary PASSCAL deployment in Turkey
right after the 1999 Mw7.4 Izmit and Mw7.1 Duzce earthquakes ( Ben-Zion et al., GJI, 2003) , and play with them. The seismograms shown below were also used in Peng & Ben-Zion (GJI, 2004) to document crustal anisotropy near the Karadere-Duzce branch of the North Anatolian Fault that ruptured during those two earthquakes.
% sac
SAC> r YJ.BV.EHZ.SAC # read the vertical-component data recorded by station BV`
SAC> p # plot the data
SAC> r YJ.BV.EH?.SAC # read all three component data recorded at station BV
SAC> p1 # plot multiple data in one page
SAC> xlim 0 7 # zoom in at 0 - 7 s relative to the waveform reference time (typically the origin time)
SAC> p1 # you have to plot it again to reflect the change
SAC> xlim a -2 t0 +4 # zoom in at 2 s before P, and 4 s after S arrivals
SAC> ylim all # enforce that y axis has the same scale
SAC> p1 # you have to plot it again to reflect the change
SAC> bd sgf # Begins plotting to the SAC Graphics File device driver
SAC> p1 # no X graphic is shown. But a file named f001.sgf is created
SAC> sgftops f001.sgf BV_3c.ps # Convert the SGF file into a postscript file
SAC> quit
% gs BV_3c.ps # ghostview your postscript file
- To write SAC data back to your disk, you have the following three options. See below for examples.
% cp YJ.BV.EHZ.SAC YJ.BV.EHZ.SAC.bak # make a backup of the vertical data
%sac
SAC> r YJ.BV.EHZ.SAC # read the data file
SAC> p1 # plot it out
SAC> rmean # remove the mean
SAC> bp p 2 n 4 c 2 8 # apply a two-way 2-8 Hz band-pass-filter
SAC> w over # overwrite the data back to the disk (using the original name YJ.BV.EHZ.SAC)
SAC> w append .bp # Write the data back to the disk using a new name YJ.BV.EHZ.SAC.bp
SAC> w temp.sac # Write the data back to the disk using a new name temp.sac
SAC> quit
% rm -f temp.sac BV.BHZ.sac # clean up what just created.
% cp YJ.BV.EHZ.SAC.bak YJ.BV.EHZ.SAC # restore the original data
- To learn the syntax of each command, use "help command".
- You can use short version of a command. For example, "p" is the same as "plot", "w" is the same
as "write".
- SAC header consists of important information about the data, such as the sampling interval, start time,
length, station location, event location, components, phase arrivals, etc. There parameters will be used by
various SAC command to process the data. It is important to keep your header information complete and updated.
- To list the SAC header information, you can open the data in SAC, and use lh (listhdr) command. To list
a specific header only, use lh header_name (e.g., lh delta). To go back to the default listing (all headers),
use lh default.
% sac # launch sac again
SAC> r YJ.BV.EHZ.SAC # read the data again
SAC> lh # list the SAC header
FILE: YJ.BV.EHZ.SAC - 1
----------
NPTS = 2409 # number of data points
B = -9.992996e+00 # begin time
E = 1.408700e+01 # end time
IFTYPE = TIME SERIES FILE # file type
LEVEN = TRUE # evenly sampled time series
DELTA = 1.000000e-02 # time increment
IDEP = VELOCITY (NM/SEC) # physical unit of the data
DEPMIN = -2.073471e+04 # minimum amplitude
DEPMAX = 1.584818e+04 # maximum amplitude
DEPMEN = 5.137106e+01 # mean amplitude
OMARKER = 0 # event origin marker
AMARKER = 1.848 # first arrival (P) marker
T0MARKER = 3.192 # t0 (S) marker
KZDATE = NOV 20 (324), 1999 # reference date
KZTIME = 00:12:55.840 # reference time
IZTYPE = GMT DAY # type of reference time
KSTNM = BV # station name
CMPAZ = 0.000000e+00 # component azimuth relative to north
CMPINC = 0.000000e+00 # component "incidence angle" reletive to the vertical
STLA = 4.075520e+01 # station latitude
STLO = 3.101490e+01 # station longitude
STEL = 2.470000e+02 # station elevation
STDP = 0.000000e+00 # station depth below surface (meters)
EVLA = 4.079930e+01 # event latitude
EVLO = 3.100330e+01 # event longitude
EVDP = 8.150000e+00 # event depth
DIST = 4.994444e+00 # source receiver distance in km
AZ = 1.686886e+02 # azimuth
BAZ = 3.486961e+02 # back azimuth
GCARC = 4.492941e-02 # great circle distance
LOVROK = TRUE # TRUE if it is okay to overwrite this file on disk
USER7 = 0.000000e+00 # User defined time picks
USER8 = 0.000000e+00 # User defined time picks
NVHDR = 6 # Header version number. Current value is the integer 6.
SCALE = 1.000000e+00 # Multiplying scale factor for dependent variable [not currently used]
NORID = 0 # Origin ID (CSS 3.0)
NEVID = 0 # Event ID (CSS 3.0)
NWFID = 2 # Waveform ID (CSS 3.0)
LPSPOL = FALSE # TRUE if station components have a positive polarity (left-hand rule)
LCALDA = TRUE # TRUE if DIST, AZ, BAZ, and GCARC are to be calculated from station and event coordinates
KCMPNM = EPZ_01 # Component name
MAG = 2.310000e+00 # Event magnitude
- A detailed description of the SAC data file format, including the header, can be found at the
SAC2000 User's manual http://www.iris.edu/manuals/sac/manual.html
- Another way to list SAC header is to use a C program saclst (originally written by Professor Lupei Zhu at SLU. The binary code is in /usr/local/geophysics/bin. This program should also be included in the latest version of SAC.
% saclst
Usage: saclst header_lists f file_lists
% saclst evla evlo stla stlo cmpaz cmpinc f YJ.BV.EH?.SAC # list the SAC header evla evlo stla stlo
YJ.BV.EHE.SAC 40.7993 31.0033 40.7552 31.0149 90 90
YJ.BV.EHN.SAC 40.7993 31.0033 40.7552 31.0149 0 90
YJ.BV.EHZ.SAC 40.7993 31.0033 40.7552 31.0149 0 0
- Like Matlab and Shell script, you can write a set of SAC commands to be executed together
in a file called SAC Macro. This makes life easy if you need to perform the same type of jobs
for many events.
- Below is a very simple Macro file that reads in three component seismograms and generate a plot
with header information in the title.
% vi plt_3c.mac
* SAC MACRO file plt_3c.mac
* reads in three component seismograms and generate a plot
* with header information in the title.
* usage:
* written by Zhigang Peng (zpeng@gatech.edu), Sun Oct 8 19:59:33 EDT 2006
* last updated on Sun Sep 13 11:58:31 EDT 2009
$KEYS sta
$default sta BV
* define keys, and default value
sc rm -f *.sgf
* remove all the previous sgf files
setbb sacfile "*.$sta$.*.SAC"
* setbb sacfile "*.$sta$.*.SAC"
setbb psfile "$sta$_3c.ps"
* set three component data
qdp off
* turn off quick dirty plot (qdp)
r %sacfile%
* read the data
rmean
* remove the mean
bp p 2 n 4 c 1 10
* use two way, 4th order bandpass butterworth filter with
* with pass band of 1 - 10 Hz
setbb ds '( concatenate ' &1,kstnm ' ' &1,kzdate ' ' &1,kztime ' ) '
evaluate to r &1,dist
title '( concatenate ' %ds% ' ', dist = ' ' %r km ') '
* get the stn name, reference time, and distance from SAC header
* and put them into the title
fileid l ul
fileid t l kcmpnm
* display the component name
* fileid Controls the file id display found on most SAC plots.
* TYPE LIST hdrlist : Define a list of header fields to display in the
* fileid.
* LOCATION UL : Place file id in upper left hand corner.
title on
xlabel "Time @(s@)"
ylabel "Velocity @(nm/s@)"
* set the xlabel and ylabel
ylim all
* set the y-axis to be the same
xlim o o +8
* set the time axis to be from origin time to 8 sec after the origin time
p1
* plot on the x window first
bd sgf
* start the SAC Graphics File device driver
p1
* plot it again to f001.sgf
sgftops f001.sgf %psfile%
* convert it to postscript, done
sc rm -f *.sgf
* remove all the previous sgf files
sc gs %psfile%
* view the outout eps file
quit
* quit SAC
- To run the SAC Macro file, use one of the following three ways:
%sac
SAC> m plt_3c.mac sta BV
# or
% printf "m plt_3c.mac sta BV\nq\n" | sac
# or
% sac plt_3c.mac sta BV
- Again, a detailed description of the SAC Macro can be found online at the SAC2000 User's manual
http://www.iris.edu/manuals/sac/SAC_Manuals/SACMacros.html
- Many basic tools are available in the SAC program. Below I will show a few examples of
these tools and functions. A complete list of the functions can be found online at SAC Home Page .
-
Rotate a pair of data.
We can easily rotate a pair of data components through an angle in SAC using the "rotate" command.
Each pair must have the same station name, event name, and sampling rate.
The most common practice is seismology is to rotate two horizontal component seismograms to be along
and perpendicular to the "great circle path" (i.e., radial and transverse component). This means that
CMPAZ must be defined (0 degrees for the North component and 90 degress for the East component)
and that CMPINC must be 90 degrees. See below for an example of rotating
two horizontal seismograms recorded at the broadband station BK.PKD near the Parkfield section of the
San Andreas fault and generated by the 2002/11/03 Mw7.9 Denali Fault earthquake in Alaska. Waveforms
recorded at this and other nearby stations are used to study non-volcanic tremor in central California
triggered by teleseismic events (Peng et al., GRL, 2008; Peng et al., JGR, 2009).
Note: the rotate command requires that the two horizontal component data to be the same length. Use the cut command to cut them into the same length if needed. Also
please make sure that the sac headers CMPINC and CMPAZ have the correct values.
% sac
SAC> r BK.PKD.HHN.SAC # read the north component seismogram
SAC> ch CMPAZ 0 CMPINC 90 # change the incident and azimuth of the instrument to be north and horizontal
SAC> wh # save the header (the above three lines are not needed if the headers are already in)
SAC> r BK.PKD.HHE.SAC # read the east component seismogram
SAC> ch CMPAZ 90 CMPINC 90 # change the incident and azimuth of the instrument to be east and horizontal
SAC> wh # save the header (the above three lines are not needed if the headers are already in)
SAC> r BK.PKD.HHN.SAC BK.PKD.HHE.SAC # read two horizontal component seismograms
SAC> rotate to GCP # rotate to the "great circle path"
SAC> w BK.PKD.HHR.SAC BK.PKD.HHT.SAC # save the rotated data (first radial, second transverse)
SAC> r BK.PKD.HHZ.SAC BK.PKD.HHR.SAC BK.PKD.HHT.SAC # read three-component seismograms (vertical, radial and transverse)
SAC> qdp off # turn off "quick and dirty plot"
SAC> xlim 300 1500 # zoom in
SAC> p1 # plot them out, try to identify Love and Rayleigh waves
-
Filtering and spectral analysis. We can easily design different filters (either lowpass, highpass, or bandpass), and compute Fast Fourier Transform (FFT) using built SAC commands.
% sac
SAC> r BK.PKD.HHT.SAC # read the transverse-component data recorded at station PKD
SAC> rmean # remove the mean
SAC> bp p 2 n 4 c 2 8 # apply two-way 2-8 Hz band-pass-filter to the data
SAC> w append .bp # save the band-pass-filtered data as BK.PKD.HHT.SAC.bp
SAC> r more # read the original data again
SAC> p1 # plot the band-pass-filtered data (high-frequency signals from locally triggered
tremor) and the original data (mostly teleseismic surface waves from the Denali Fault earthquake)
SAC> cut 600 1600 # cut the data
SAC> r BK.PKD.HHT.SAC # read the data again
SAC> rmean # remove the mean
SAC> taper # apply a hanning window to the 0.05 width of the data
SAC> fft # compute FFT
SAC> keepam # keep only the amplitude information (and throw the phase information)
SAC> p1 # plot the amplitude spectrum
SAC> loglog # use the log-log plot
SAC> xlim 0.001 50 # set the x-axis plot range (Frequency)
SAC> p1 # plot again in log-log scale
-
Phase pick. Phase picking can be easily done in SAC via automatic phase picker "apk", or
manually phase picker "ppk".
% sac
SAC> cut 0 10 # cut the data
SAC> r YJ.BV.EH?.SAC # read three component data recorded at station BV
SAC> qdp off # turn off quick dirty plot
SAC> p1 # plot the data
SAC> ch a -12345 t0 -12345 # remove the original P and S arrival pick
SAC> apk # use the auto picker
SAC> p1 # see the auto picker result
SAC> ppk # use the manual picker
# cheat sheet,
zoom in: type "x" to define left side time window,
followed by a left click of mouse to define the right side time window
zoom out: type "o"
p arrival: type "a", or "p" at the time
where you think the p wave arrival in.
s arrival: type "t0", or "s" at the time
where you think the s wave arrival in.
quit: type "q"
SAC> ppk p 3 m on # pick all 3 component together with the same time
# to save the picked arrival, you can use the command "wh". But this command could
# mess up the data if you have used the cut command. So be careful.
SAC> wh
# save your picks, sometimes it would be useful to add "bell off" following ppk to turn off the annoying sounds!
- IRIS stands for Incorporated Research Institutions for Seismology. It is
an university research consortium dedicated to exploring the Earth's interior through the
collection and distribution of seismographic data (Ref:
http://www.iris.edu/about/).
- Basically IRIS is where seismologists get instruments to do seismic experiment,
and download seismic data for scientific research.
- Other earthquake data centers include: Southern California Earthquake Data Center
http://www.data.scec.org , Northern California Earthquake Data Center
http://www.ncedc.org ,
HiNet, Japan ,
F-Net, Japan ,
K-Net, Japan ,
KiK-Net, Japan ,
- There are many tools that are available to download seismic data from these data center. IRIS has a
nice tutorial on this: IRIS Data Management Center Data Access Tutorial .
- First, please go to the following working directory:
% cd SAC_Tutorial_WF/Denali_2002
- In this exercise, we will use WILBER II to download waveform data generated by the 2002 Mw7.9 Denali Fault earthquake.
- You can download waveforms by clicking on the following web page: http://www.iris.edu/wilber .
- Next choose Q4 2002, and the Alaska. Then you will get a page with all events with 5 deg (or other specfic range)
of the clicked point. Click on the following M6.3 event:
DATE TIME SOURCE MAG LAT LON DEPTH DESCRIPTION
2002/11/03 22:12:41.0 FARM 7.9 63.52 -147.44 4.90 CENTRAL ALASKA
- Next, you need to define the network, channel, distance, and other parameters to select the data
you would like to request from IRIS. To make the request relatively modest, let's use the following parameters:
NETWORK - AK, AZ, BK, CI, and UW, and click on PROCEED. Next, choose
CHANNEL - BHZ, BHN, BHE, DISTANCE from 0 to 180 deg.
- Next, we need to select data format, time window and user identifition. We will choose the SEED
format, 5 minutes before, and 60 minutes (1 hour) after P.
- After a few minutes, you should get an email notice saying that you data is read to pick up somewhere
in an FTP site. Please download the file and put it inside Sac_Tutorial/Denali_2002 directory.
- For the GT user, you can use the SEED data I have already downloaded from the IRIS DMC (to save time). I did not
put the SEED data online for the outside users.
- To extract SAC waveform data from the SEED volume, use the following command:
% cd SAC_Tutorial_WF/Denali_2002
% rdseed.linux -d -o 1 -f 20021102_Denali.seed
- Seismological research has changed significantly in the past few decades, partially due to the digital
revolution. In most case, we will not be able to publish a paper based on only one or two seismograms (I was
told that people in the early 1950-1960s can do that). In addition, some of the most exciting findings within
seismology come from the analysis of massive continous seismic waveforms that are
otherwise considered as noise (e.g., non-volcanic tremor, ambient noise tomography, etc).
Hence, we need to analyze hundreds and thousands
(sometimes even more) of seismic data first. How can we do that?
- Basically we have to reply on computers to do most of the dirty job. 8-) Shell script language (e.g., gawk,
sed, Makefile and Perl) will help you a lot in this case. Below are some examples to show
how you can interact with SAC using shell script.
- Example 1: rename the SAC filename to something like NET.STN.BH[ZNE].SAC
% mkdir backup
% ls *.SAC | gawk -F"." '{print $7,$8}' | sort | uniq -c | gawk '{if ($1 == 3) print $2,$3}' | sort | uniq > stn_3c.id
# find the station name which has three seismograms (three components)
% ls `gawk '{print "*"$1"."$2".*SAC"}' stn_3c.id` | gawk -F"." '{print "mv "$0,$7"."$8"."$10"."$12}' | sh
# rename the SAC filename to NET.STN.BH[ZNE].SAC
% mv 2002*.SAC backup
# move the rest into backup for now
- Example 2: put earthquake and station component information into the SAC header.
% gsact 2002 11 03 22 12 41 518 f *.SAC | gawk '{print "r "$1"\nch o "$2"\nwh"} END{print "quit"}' | sac
# put the eq. origin time into the SAC header
% printf "r *.SAC\nch evlo -147.45290 evla 63.51410 evdp 4.2 mag 7.9\nwh\nq\n" | sac
# put the earthquake location information into the SAC header
% saclst o f *.SAC | gawk '{printf "r %s\nch allt %.5f\nwh\n",$1,$2*(-1)} END{print "quit"}' | sac
% shift the origin time to be at 0 s
% printf "r *.BHE.SAC\nch cmpaz 90 cmpinc 90\nwh\nq\n" | sac
% printf "r *.BHN.SAC\nch cmpaz 0 cmpinc 90\nwh\n\q\n" | sac
% printf "r *.BHZ.SAC\nch cmpaz 0 cmpinc 0\nwh\n\q\n" | sac
# put the component information into the SAC header
- Example 3: cut the two horizontal component data so that they will have the same length, which are
needed for rotation in SAC.
% saclst e f `gawk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` # list the end of the two horizontal component data
% saclst e f `gawk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` |\
paste - - | gawk '{if ($2> $4) print $4;else print $2}' > BH_EN_end_time.dat
% saclst b f `gawk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` |\
paste - - | gawk '{if ($2< $4) print $4;else print $2}' > BH_EN_start_time.dat
# find out the minimum of the end time and save it
% paste -d" " stn_3c.id BH_EN_start_time.dat BH_EN_end_time.dat |\
gawk '{print "cut "$3,$4"\nr "$1"."$2".BH[EN].SAC\nrmean\nrtrend\nrotate to GCP\nw "$1"."$2".BHR.SAC "$1"."$2".BHT.SAC"} END {print "q"}' | sac
# cut the two horizontal component data, remove mean and trend,
# rotate them to great circle path, save the output
% gawk '{print "r "$1"."$2".BH[ZNE].SAC\nrmean\nrtrend\nw over"} END {print "q"}' stn_3c.id | sac
# remove mean and trend for the original three component data.
- An important step in data processing is to document what we have done to the data, so that other
people can reproduce your results. Also sometimes you need to reprocess the data, but maybe you have
forgot the commands and parameters you have typed a few months ago.
- There are two ways to document your scripts. The first way is simply to write a shell script.
This could be assigned as a potential homework for you. You can define different steps using flags,
and simply copy the command line from the above examples into each step.
- Another way is to write a Makefile, and put each step seperately. Makefile was originally used for compiling
a collection of computer codes. It can be used to store simple commands if puting them in a seperate shell script
file is not needed.
- See Makefile in your working directory. For example, if you want to redo step 3, simply type "make step3"
in your command line, it will automatically execute what is written underneath the flag step3. If you want to clean
up all you have done here, simply type "make clean". If you want to reprocess all the steps, type "make all".
- Note that to make the Makefile work, you need to change "$" to "$$" when you refer to variables in gawk.
- Additional tutorial on Makefile can be found online, such as Make a tutorial .
- To convert between SAC binary and ASCII files, the easiest way is to use SAC's built command 'write alpha' to save the data into ASCII format
SAC> r AZ.RDM.BHZ.SAC
SAC> w alpha AZ_RDM_BHZ.txt
- An alternative would be to use the following command lines developed by us (originally from Dr. Lupei Zhu at SLU.
% sacdump AZ.RDM.BHZ.SAC | head
% sacdump AZ.RDM.BHZ.SAC > AZ_RDM_BHZ.dat
% # the output is two ASCII columns with the first being the time and second being the actual data.
% sacdump_slice 0 2000 AZ.RDM.BHZ.SAC > AZ_RDM_BHZ_cut.dat
% # in this code, you can pre-define the time window where you would like the data to be converted into the two column ASCII file.
% gawk '{print $2}' AZ_RDM_BHZ_cut.dat | col2sac AZ.RDM.BHZ.cut.SAC 0.025 0
% # the command col2sac convert the one column data from stdin into a sac file, usage: col2sac sac_file_name delta_t t0
- We can plot seismic data in SAC, but the output figure does not look great. If you are a GMT fan like me,
you may wonder if there is a way we can plot SAC data using GMT. The answer is yes.
- Dr. Lupei Zhu at SLU wrote a code called "pssac"
that can do the job. I already put the binary code in /usr/local/geophysics/bin. The source code of "pssac" can be
downloaded from the following link: GMT-style SAC trace plotting pssac and the new one for GMT4.5.
- The syntax of pssac is very similar to a typical GMT command. Below are two simple examples.
The first one plots a record section of vertical component seismograms for the 2006/10/15 Hawaii earthquake.
The second plot the three component and radial and transverse component data at station GSC.
% gmt4 # command line to use the GMT version 4 (which can run pssac)
% pssac -R0/2000/5/38 -JX6i/9i -Edt-3 -M0.3 -K -P -Ba200f100:"Time (s)":/a5f1":Distance (deg)":WeSn *.BHZ.SAC > Denali_BHZ_record_section.eps
% printf "0 2000\n" | gawk '{printf "%f %f\n%f %f\n",$1,$1*4,$2,$2*4/111.19}' | psxy -R -JX -K -O -P -W1tap/255/0/0 >> Denali_BHZ_record_section.eps
% printf "1000 39 20 0 4 MC 11/02/2002 Denali Fault Earthquake\n" | pstext -R -JX -O -P -N -G0/0/0 >> Denali_BHZ_record_section.eps
% gs Denali_BHZ_record_section.eps
% pssac -R300/1800/0/6 -JX6i/4i -Ent-3 -M1 -K -P -Ba500f100:"Time (s)":S AZ.RDM.*.SAC > Denali_AZ_RDM_BH.eps
% ls AZ.RDM.*.SAC | gawk -F"." '{printf "%f %f 20 0 4 ML %s %s\n",300,NR+0.2,$2,$3}' | pstext -JX -R -O -P >> Denali_AZ_RDM_BH.eps
% gs Denali_AZ_RDM_BH.eps
- You can also type make Denali_BHZ_record_section.eps Denali_AZ_RDM_BH.eps to generate these two plots.
- If you like to use Matlab, you may want play with SAC data in Matlab. I have developed a set of
matlab subroutines named MatSAC that will read and write SAC data in Matlab. These subroutines are in MatSAC
directory. To read a SAC file into Matlab, use the following syntax:
[t,data,SAChdr] = fget_sac(filename)
- The MatSAC package can be downloaded from the following website: MatSAC .
After that, it is best to put it under ${home}/matlab directory.
cd ~/matlab
wget http://geophysics.eas.gatech.edu/people/zpeng/Teaching/SAC_Tutorial/MatSAC.tar.gz
tar zxvf MatSAC.tar.gz
vi startup.m
path(path,'~/matlab/MatSAC');
- There is a simple matlab script named gen_bp_PKD.m in your working directory.
It will use fget_sac to read transverse component data at station PKD,
and plot both the broadband and band-pass-filtered data.
cd SAC_Tutorial/SAC_Tutorial_WF # go to the first tutorial example
matlab # launch matlab, with the fancy interface, or better use the command line only
matlab -nojvm # launch matlab without the java machine
>> gen_bp_PKD # run the matlab script, which calls fget_sac command.
Part 3: Advanced SAC tools
- Please use the following examples that are recorded by station NC.CCO and generated
by a repeating cluster along the Calaveras fault
in central California. This station recorded clear fault zone head waves that refract along the bi-material
fault interface, and has been studied in Zhao and Peng (GRL, 2008) . We can perform waveform cross-correlation in SAC using the following
commands:
cd SAC_Tutorial/SAC_Tutorial_WF
SAC> cut 1 2 # cut the data at 1-2 s after the origin time (around the P waves)
SAC> r *.NC.CCO.EHZ.SAC # read all the data
SAC> rmean # remove the mean
SAC> taper # apply a cosine taper with 5% lenth
SAC> correlate # do the cross-correlation with the first trace
SAC> xlim -0.2 0.2 # zoom in around the zero lags
SAC> p1 # plot them to see the time shifts
SAC> p2
- Another example is to use the sac_wfcc command (modified from the code written by Dr. Lupei Zhu at SLU)
to align the P wave (or any other phases), and compute the correlation coefficients. For outside users, the
source code for sac_wfcc can be downloaded at the following link sac_wfcc.tar.gz .
% saclst a f *.NC.CCO.EHZ.SAC | sac_wfcc -D-0.2/0.8/0.5
# list the original P wave arrival, and perform waveform cross-correlations using 0.2 s before
# and 0.8 s after the P arrival. The output is the wf, the P arrival after time shift to best align
# with the first trace, and the correlation coefficient.
% saclst a f *.NC.CCO.EHZ.SAC | sac_wfcc -D-0.2/0.8/0.5 | gawk '{print "r "$1"\nch t6 "$2"\nwh"} END {print "q"}' | sac
# put the aligned P wave arrival into the SAC header t6
% saclst t6 f *.NC.CCO.EHZ.SAC | gawk '{printf "r %s\nch allt %.5f\nwh\n",$1,$2*(-1)+5} END {print "q"}' | sac
# shift the trace so that they are aligned at 5 s.
SAC> cut 4.5 5.5
SAC> r *.NC.CCO.EHZ.SAC
SAC> qdp off
SAC > p1
SAC> p2
# as you see, the waveforms after alignment are more coherent than the original data.
- Stacking helps to enhance the coherent signal and suppress the non-coherent random noise. Many
advanced array processing tools are based on the idea of stacking. Here we only use the SAC command
to perform direct stacking (i.e., no time shifts). There are many existing tools available for advanced
array analysis (e.g., stant stack).
SAC> cut 4 6 # cut the data from 4 to 6 s (around the P waves)
SAC> r *.NC.CCO.EHZ.SAC # read all the data
SAC> rmean # remove the mean
SAC> w append .P # save the cut data
% ls *.NC.CCO.EHZ.SAC.P | gawk '{if (NR == 1) print "r "$1;else print "addf "$1} END {print "div "NR"\nw CCO.stack.SAC.P\nq"}' | sac
# read the first data, add the rest data, divide with the total number, and save the output
SAC> r *.NC.CCO.EHZ.SAC.P # read the P wave data
SAC> color black inc # set the color to start from black
SAC> p1 # plot them
SAC> p2
Under development. To be updated in Spring 2013.
|