Introduction to Seismic Analysis Code (SAC)

by Zhigang Peng

This website contains a brief tutorial on Seismic Analysis Code (SAC). It is part of a mini course offered to the Geophysics graduate students at GT in Fall 2006 and 2009 by Professor Andy Newman and me. The tutorial consists of the following sections:

SAC overview

How to run it at GT geophysics laboratory

  • SAC is already installed at "/usr/local/geophysics/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/aux
% echo $PATH
  • You should have SACAUX defined as "/usr/local/geophysics/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/aux
export PATH=/usr/local/geophysics/bin:$PATH

SAC data format

  • The standard data format for SAC is binary. A binary SAC file contains a fixed length ASCII data that describe the following binary data.
  • For seismic data this means a single data component recorded at a single seismic station. SAC does not currently work on multiplexed data. The data will generally be evenly spaced time series 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), while in LINUX (Intel processors), it is stored according to the "little-endian" byte order.
  • Since all of us are working in LINUX now, make sure you convert the data first 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.

Read/Write/Plot seismic data in SAC

  • Let's start with an example. 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/Sac_Tutorial_2006.tar.gz 
% cd Sac_Tutorial_2006/SAC_WF
  • Next, we will read some seismograms, and play with it. 
    % sac
SAC> r BV.z		# read the data file
SAC> p			# plot a single data
SAC> r BV.?		# 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 start 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. 
    %sac
SAC> r BV.z		# read the data file
SAC> p1			# plot it out
SAC> rmean		# remove the mean
SAC> w over		# overwrite the data back to the disk (using the original name BV.z)
SAC> w append .sac	# Write the data back to the disk using a new name BV.z.sac
SAC> w temp.sac		# Write the data back to the disk using a new name temp.sac
SAC> quit
% rm -rf temp.sac BV.z.sac # clean up what just created.
  • 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

  • SAC header consists of important information about the data, such as sampling interval, start time, length, station location, event location, components, phase arrival, 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 BV.z			# read the data again
SAC> lh				# list the SAC header
  FILE: BV.z - 1
 ----------

    NPTS = 2409			# number of data points
       B = -9.676000e+00	# begin time
       E = 1.440400e+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.317		# event origin marker
 AMARKER = 2.165		# first arrival (P) marker
T0MARKER = 3.509		# t0 (S) marker
  KZDATE = NOV 20 (324), 1999	# reference date
  KZTIME = 00:12:55.523		# 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.llnl.gov/sac/SAC_Manuals/UserManualOverview.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. 
    % saclst
Usage: saclst header_lists f file_lists
% saclst evla evlo stla stlo f BV.?				# list the SAC header evla evlo stla stlo
BV.e      40.7993     31.0033     40.7552     31.0149
BV.n      40.7993     31.0033     40.7552     31.0149
BV.z      40.7993     31.0033     40.7552     31.0149

SAC Macro file

  • 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, Sun Oct  8 19:59:33 EDT 2006

$KEYS sta
$default sta BV

sc rm -f *.sgf
* remove all the previous sgf files

* define keys, and default value

* setbb sacfile "$sta$.[zne]"
setbb sacfile "$sta$.*"

setbb psfile "$sta$_3c.eps"
* set three component data

qdp off
* turn off quick dirty plot 8-)

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\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.llnl.gov/sac/SAC_Manuals/UserManualOverview.html

Other basic tools available in SAC

  • 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 Commands Listing: Functional .

  • 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 and that CMPINC must be 90 degrees. See below for an example of rotating two horizontal seismograms at Kota Tinggi in Malaysian and generated by the 2006/07/17 M7.7 JAVA Tsunami earthquake. 
    % sac
SAC> r KOM.BHN.SAC KOM.BHE.SAC			# read two horizontal component seismograms
SAC> rotate to GCP				# rotate to the "great circle path"
SAC> w KOM.BHR.SAC KOM.BHT.SAC			# save the rotated data (first radial, second transverse)
SAC> r KOM.BHZ.SAC KOM.BHR.SAC KOM.BHT.SAC	# read three-component seismograms
SAC> qdp off					# turn off "quick and dirty plot"
SAC> xlim 100 1000				# zoom in 
SAC> p1						# plot them out, try to identify Love and Rayleigh waves
SAC> cut 550 650				# prepare the data window for Rayleigh wave
SAC> r KOM.BHZ.SAC KOM.BHR.SAC			# read the vertical and radial component only
SAC> xlim off					# turn off the zoom in time window
SAC> p2						# plot them on top of each other
SAC> ppm					# plot the partical motion
SAC> quit					# done

  • Phase pick.

    Phase picking can be easily done in SAC via automatic phase picker "apk", or manually phase picker "ppk". 
    % sac
SAC> r BV.?			# read three component data
SAC> qdp off			# turn off quick dirty plot
SAC> p1				# plot the data
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
SAC> wh				# save the header information, be careful

  • Compute cross correlation.

    Below is an example of how to compute cross correlation for waveforms generated a group of repeating earthquakes and recorded at station CCO in northern California. 
    %sac
SAC> cut t0 -0.5 t0 +2.5	# cut the data 0.5 s before and 2.5 s after the s arrival
SAC> r CCO*.EHZ			# read all the waveforms
SAC> p2				# plot the data on top of each other
SAC> correlate			# compute waveform cross correlation with respect to the first data
SAC> xlim -0.5 0.5		# zoom in the cross-correlation functions
SAC> p1				# plot out the cross-correlation functions

Download seismic data from IRIS (and other data center)

  • 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 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/Sac_Tutorial_2006.tar.gz 
% cd Sac_Tutorial_2006/HAWAII_20061015
  • In this tutorial, we will use WILBER II to download waveform data generated by the recent Hawaii M6.3 earthquake. There is a report on this event at IRIS Special Event Pages, where there is link to view data via Wilber II.
  • But let's follow the standard procedure, and click on the following web page: http://www.iris.edu/cgi-bin/wilberII_page1.pl .
  • Next click on Hawaii Island. 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                    
2006/10/15 17:07:48.3 SPYDER® 6.3  19.84 -156.06  24.00 HAWAII
  • 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 - ALL, CHANNEL - BHZ, BHN, BHE, DISTANCE from 0 to 180 deg, and check DISTRIBUTE about every 1 deg (so we will not get too many data for this tutorial). You can click on each station to preview waveforms.
  • Next, we need to select data format, time window and user identifition. We will choose the default SEED format, 5 minutes before, and 240 minutes (4 hours) 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_2006/HAWAII_20061015 directory.
  • In case there is a problem, you can download the data I have requested online at ftp://ftp.iris.washington.edu/pub/userdata/Zhigang_Peng/HAWAII_20061015_M63. This file will only be there for a few days, though. Or you can copy it from /usr/local/geophysics/sac/HAWAII_20061015_M63.seed
  • To extract SAC waveform data from the SEED volume, use the following command: 
    % rdseed -d -o 1 -f HAWAII_20061015_M63.seed

Use shell script (e.g., gawk, sed, Makefile, Perl) to interact with SAC

  • Seismological research has changed significantly in the past few decades, partially due to the digital revolution. In most case, you 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). Instead, you 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. 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 STN.BH[ZNE].SAC 
    % mkdir backup
% saclst b e f 2006*.SAC | gawk '{if ($3-$2<1000) print "mv "$1" backup"}' | sh
# move those sac files which has data length less than 1000 s
% ls *.SAC | gawk -F"." '{print $8,$9}' | sort | uniq -c | gawk '{if ($1 == 3) print $2}' | sort | uniq > stn_3c.id
# find the station name which has three seismograms (three components)
% ls `gawk '{print "*."$1".*SAC"}' stn_3c.id` | gawk -F"." '{print "mv "$0,$8"."$10"."$12}' | sh
# rename the SAC filename to STN.BH[ZNE].SAC
% mv 2006*.SAC backup #move the rest SAC file to backups
  • Example 2: 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".BH[EN].SAC"}' stn_3c.id` # list the end of the two horizontal component data
% saclst e f `gawk '{print $1".BH[EN].SAC"}' stn_3c.id` |\
% paste - - | gawk '{if ($2> $4) print $4;else print $2}' > BH_EN_end_time.dat
# find out the minimum of the end time and save it
% paste -d" " stn_3c.id BH_EN_end_time.dat |\
gawk '{print "cut 0",$2"\nr "$1".BH[EN].SAC\nrmean\nrtrend\nrotate to GCP\nw "$1".BHR.SAC "$1".BHT.SAC"}' | sac
# cut the two horizontal component data, remove mean and trend,
# rotate them to great circle path, save the output
% gawk '{print "r "$1".BH[ZNE].SAC\nrmean\nrtrend\nw over"}' 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. I wrote a simple shell script named preprocess.bash in your working directory that will do the whole job. Note that sometimes you may want to redo only part of the job. So I add a flag for each step, and you can reset the flag to turn that step on and off.
  • 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's 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 "$$".
  • Additional tutorial on Makefile can be found online, such as Make a tutorial .

Plot SAC data in GMT

  • 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 .
  • 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. 
    % pssac -R-100/5000/-5/180 -JX6i/9i -Edt-3 -M0.3 -K -P -Ba1000f100:"Time (s)":/a50f10":Distance (deg)":WeSn *.BHZ.SAC > HAWAII_BHZ_record_section.eps
% printf "2550 170 20 0 4 MC 10/15/2006 HAWAII Earthquake\n" | pstext -R -JX -O -P -G255/0/0 >> HAWAII_BHZ_record_section.eps
% gs HAWAII_BHZ_record_section.eps
% pssac -R100/2000/0/6 -JX6i/6i -Ent-3 -M1 -K -P -Ba1000f100:"Time (s)":S *GSC*.SAC > HAWAII_GSC_BH.eps
% ls GSC.*.SAC | gawk -F"." '{printf "%f %f 20 0 4 ML %s\n",300,NR+0.2,$2}' | pstext -JX -R -O -P >> HAWAII_GSC_BH.eps
% gs HAWAII_GSC_BH.eps
  • You can also type make HAWAII_BHZ_record_section.eps HAWAII_GSC_BH.eps to generate these two plots.

Read/Write/Plot SAC data in Matlab

  • 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
  • There is a simple matlab script named plot_fft_BHZ.m in your working directory. It will use fget_sac to read vertical component data at station GSC, zoom in around P and Rayleigh waves, compute FFT, and plot the amplitude spectra. Please run it in Matlab, and Feel free to modify it as needed.

Future development

  • Removing instrument response using SAC
  • Performing waveform cross-correlation using SAC
Last updated by Zhigang Peng Wed Aug 20 09:12:26 EDT 2008