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

SAC overview


How to run it at GT geophysics laboratory

  • 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)
    

SAC data format

  • 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.

Read/Write/Plot seismic data in SAC

  • 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

  • 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
    

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 (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

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 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

Various misc codes that can be used to 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

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 go to the following working directory:
    % cd 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 .
  • WILBER II is scheduled to be retired on September 30, 2013. It will be replaced by Wilber 3 . I will update the tutorial later this year.
  • 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.
  • To save time, you can use the SEED data I have already downloaded from the IRIS DMC.
  • To extract SAC waveform data from the SEED volume, use the following command:
    % cd Denali_2002
    % rdseed -d -o 1 -f 20021103_Denali.seed
    

Use shell script (e.g., awk, 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, 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 .

Convert data between SAC binary and ASCII files

  • 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
    

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 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.

Plot SAC data in record section in SAC

  • 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

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 . 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

Waveform cross-correlations

  • 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

  • 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
    

    Removing instrument responses

    Under development. To be updated in Spring 2014.


Last updated by Zhigang Peng Mon Aug 5 23:21:10 EDT 2013