Introduction to Seismic Analysis Code (SAC), 2013 beta version
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:
Part 1: SAC Basics
- 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_2013.tar.gz
% cd SAC_Tutorial/SAC_Tutorial_WF
For USArray short course user, please untar the following file to your home directory,
or any subdirectory where you would like to put the tutorial material.
% tar zxvf /Volumes/USArray/presentations/SAC_Tutorial_2013/Sac_Tutorial_2013.tar.gz
% cd SAC_Tutorial/SAC_Tutorial_WF
% # or if you would like to copy the tar.gz to your working directory first, and then untar.
% cp /Volumes/USArray/presentations/SAC_Tutorial_2013/Sac_Tutorial_2013.tar.gz .
% tar zxvf Sac_Tutorial_2013.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_2013.tar.gz
% tar zxvf Sac_Tutorial_2013.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
% gs -c '-100 600 translate 270 rotate' -f BV_3c.ps # ghostview by rotating it into landscape, may not work
- 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
SAC Users Guide v101.6 - June 2013, and SAC Data Format .
- 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. Outside users can download the source code at
sac_msc.tar.gz" . 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"
setbb pdffile "$sta$_3c.pdf"
* 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
sc 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
sc ps2pdf %psfile% %pdffile%
* convert the ps to pdf 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 bell off # pick all 3 component together with the same time, turn the annoying bell off
SAC> wh
# 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.
Part 2: Interact with SAC
- Below please find the link to download various source codes (mostly in C) that have been used by us (and many others) to
interact with SAC. Many of the codes were originally developed by scientists/students at Caltech, in particular, Dr. Lupei Zhu,
Professor of Geophysics at St. Louis Univ. Some of these programs are used in this tutorial. See below for a list of the codes
available (and Disclaimer) in the package.
- sac_msc readme
- sac_msc.tar.gz
- 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., awk,
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 | awk -F"." '{print $7,$8}' | sort | uniq -c | awk '{if ($1 == 3) print $2,$3}' | sort | uniq > stn_3c.id
# find the station name which has three seismograms (three components)
% ls `awk '{print "*"$1"."$2".*SAC"}' stn_3c.id` | awk -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 | awk '{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 | awk '{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 `awk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` # list the end of the two horizontal component data
% saclst e f `awk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` |\
paste - - | awk '{if ($2> $4) print $4;else print $2}' > BH_EN_end_time.dat
% saclst b f `awk '{print $1"."$2".BH[EN].SAC"}' stn_3c.id` |\
paste - - | awk '{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 |\
awk '{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
% awk '{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 awk.
- 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 commands (sacdump and sacdump_slice) 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.
% awk '{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.
- Please note that right now pssac is not compatible with the latest version of GMT. Dr. Zhu and I will contact Dr. Paul Wessel who manages
GMT to see if there is a way to include pssac in the new release of GMT. Please stay tuned!
- 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, only for GT user), or type 'source /usr/local/geophysics/etc/gmt4rc'
% 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" | awk '{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 | awk -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 the GMT pssac command did not work, you can also plot SAC data in record section using SAC's Signal Stacking Subprocess (SSS) package . Below is an example.
SAC> r CI.*.BHZ.SAC # read the vertical component data recorded by the CI network in California
SAC> sss # enter the the Signal Stacking Subprocess
SAC> plotrecordsection # plot the record section (note that the corresponding event and station headers are needed
SAC> timewindow 0 3000 # change the time window
SAC> vm 1 refractedwave vapp 3.5 # A reduction velocity of 3.5 km/s (Rayleigh waves) can be plotted
SAC> quitsub # quit the SSS package
- See more tips/examples at Peter Shearer's website http://http://mahi.ucsd.edu/shearer/227C/sac.txt
- 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');
# or download the zip file at http://geophysics.eas.gatech.edu/people/zpeng/Teaching/SAC_Tutorial/MatSAC.zip
- 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.
- Below please find a link to another matlab script called gen_spectrogram.m . It will generate the third panel that shows the spectrogram (i.e., frequency versus time) of the recorded signal (mainly teleseismic P wave, its coda, and triggered tremor at Parkfield). It needs a subroutine called eqfiltfilt.m , which is a two-pass, zero-phase butterworth filter.
cd SAC_Tutorial/SAC_Tutorial_WF # go to the first tutorial example
wget http://geophysics.eas.gatech.edu/people/zpeng/Software/matlab/gen_spectrogram.m
wget http://geophysics.eas.gatech.edu/people/zpeng/Software/matlab/eqfiltfilt.m
matlab -nojvm # launch matlab without the java machine
>> gen_spectrogram # run the matlab script to compute spectrogram, which calls fget_sac command.
>> gen_spectrogram('PKD','HH','T','BK'); # or be specific about the parameters
- Finally please find a link to another matlab script called sac2wav.m . It will read the raw seismic data, apply a 0.5 Hz high-pass filter (See Peng et al. (SRL, 2011) for why doing this), and convert the seismic data into audible sounds. If it works, you will generate an audible
seismogram like this BK.PKD.HHT.SAC.wav .
cd SAC_Tutorial/SAC_Tutorial_WF # go to the first tutorial example
wget http://geophysics.eas.gatech.edu/people/zpeng/Software/matlab/sac2wav.m
matlab -nojvm # launch matlab without the java machine
>> sac2wav # run the matlab script to convert SAC to wave, which calls fget_sac command.
>> sac2wav('BK.PKD.HHZ.SAC',100,1,0,2000); # or be specific about the parameters
- For more information, please see Earthquake Music by Zhigang Peng .
- In addition, IRIS recently made a data product called SeisSound, the Audio/Video Seismic Waveform Visualization . It has a SeisSound Repository , which contains various animations with sounds made from USArray and other openly available data from IRIS, and a source code if you would like to make them by yourselves.
- Finally, IRIS DMC allows you to download the audio wav format of any seismic data using webservices. Please click the following link for URL Builder: timeseries and choose the output to be 'audio'. The click the link on that page, or below to download the audio files.
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 | awk '{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 | awk '{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 | awk '{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 2014.
|