Some SAC Tips & Macro Examples









CJA's Guide Index | Ammon's Home | Department of Geosciences


What is SAC?

SAC is the the Seismic Analysis Code developed at Lawrence Livermore National Laboratory over the last 10-15 years. The package is widely used in the seismological community and many data are available in SAC format. See the manuals that are distributed around the department for details. In this page, I simply collected a few hints at getting the most out of the package, and details that are always useful, but often forgotten (such as the syntax for the library calls to read and write SAC files).

Another useful place for information is the online SAC information at Livermore's SAC Web Pages. They have an online manual and other useful information at that site. These notes are no substitude for the manual - they are designed to help your find some syntax or see some more complex macro examples.

Contents


SAC Data

The SAC data format is binary and can be written using the SAC library which is a part of the SAC package. Briefly, each SAC file consists of a header that stores information about the time, sampling, etc. of the signal and information on the location of the station and event. Each file contains one seismogram which makes SAC a little unwieldy for large array studies, but fine for many seismic applications.

Although the SAC data format is very rigid and in binary format, you can read ascii data directly into SAC using the readalpha [ra] command. For example, you can read two columns of ASCII data using the command




sac> ra content p input_file_name

Or, if the first column was the x-value, and the next three columns were the y-values of three different signals, use:




sac> ra content xy3 input_file_name

You need to set the header values for these data, but the approach is great for quickly examining numerical data using SAC's plotting capabilities.

Contents


Reading and Writing SAC files from Your Programs

The SAC library contains a suite of routines for writing and reading SAC files from you own computer programs. You can call these programs if you link your programs with the SAC library. See the SAC manual for details, these are just some statements for a quick reminder.

The simplest way to read an evenly spaced SAC file is using the routine rsac1:




call rsac1(file_name,seis,npts,begin_time,dt,maxpts_allowed,error_flag) 
if(error_flag .gt. 0) then
   write(stdout,*)'Fatal error reading:', file_name
   stop
else if(error_flag .lt. 0) then
   write(stdout,*)'Warning - problem reading:', file_name
end if

The simplest way to write a SAC file is using wsac1:




call wsac1(file_name,seis,npts,begin_time,dt,error_flag)
if(error_flag .ne. 0) then
   write(stdout,*)'Problem writing:', file_name
end if 

You can also set and access all the header variables in the SAC file if you use the routines rsac0 and wsac0 instead of rsac1 and wsac1 - see the manual.

Contents


SAC Plot Files

To generate a plot in SAC you must use the SAC program. The program generates files known as SAC plot files (f???.sgf). As you run the program, the number (???) increases from 001 to 999 (if you generate that many). You can control the specific names of the ouput files (see the manual). Each time you start SAC, the graphics file counter is set to 001, so you will overwrite plots if they exist in the directory. IF a graphics file is important, rename the file to something more descriptive.

The f???.sgf files are binary files and must be processed using one of several programs supplied with the SAC program. The two most common are the sgfplot and sgftops. The first plots a file on the computer screen




unix> sgfplot f001.sgf x 

will display the file under X-windows. The second program, sgftops, converts the sgf commands into postscript, which is the language that the laser printers use to render graphics and text. To plot the sgf files, use a command such as:




unix> sgftops f001.sgf my_plot.ps

Then you can plot the postscript file using




unix> lpr -Plj my_plot.ps

Or you can transfer the file into Adobe Illustrator for editing, adding a figure caption, etc. on the Macintosh (See the notes on GMT on how to make that transfer). The plot contains an identification string in the lower left that lists the date and the file path, which is very useful for figureing out where a figure came from...

You can vary the line thickness options using sgftops as well. The complete usage information for the program is:




usage: sgftops sgf_file ps_file [ line_width id_flag scale_flag]
 where:  line_width = 1, 2, or 3
         id_flag    = yes or no
         scale_flag = yes or no
  example: sgftops foo.sgf foo.ps 2 yes no
      for line width 2, id label, no scaling
Contents


SAC Macros

Without a doubt, the feature of SAC that makes it very useful is the macro programming ability. See the front of the manual for the syntax in SAC macros. However, some simple examples are presented below for your reference. In essence, SAC macros execute SAC commands, but they can also be used for calculations involving information stored in the SAC header and for automatically processing entire directories full of files.

Contents


Some quick hints for SAC macros:

To set up a blackboard or temporary variable:

setbb my_variable 25

To access a header value of the first file

&1,depmax

For the third file:

&3, depmax

To perform a calculation:

evaluate to my_answer &1,depmax * 1.1

To use a value of a Blackboard variable (call the s1 ):

ch b %s1

To construct strings, append a $ on each end of a string argument (use the % sign for balackboard strings)

r $1$_dsp.z

Three example macros follow. They are not the most efficient set of macros, but they were relatively easy to write (much easier than writing a program).

Contents


Example 1: Generating a Simple Plot

Here is a very simple macro that reads in three seismograms (vertical, radial, and transverse) and makes a plot with the date, distance and azimuth in the title. First, I set up some graphics parameters and then I read in the vertical component to get the header information. I use this information and some balcboard variables to compose the plot title. Then I read in the three components, apply a Butterworth filter, and generate the plot. Here is an image of a figure generated using this macro: example output.





xvp .2 .8
fileid l ul
fileid t l kcmpnm
title on
rh $1$.z
*
evaluate to c &1,dist + 0.5
evaluate to r ( INT %c )
setbb rs ( before '.' ' %r ' )
evaluate to c &1,az + 0.5
evaluate to a1 ( INT %c )
setbb a2 ( before '.' ' %a1 ' )
setbb ds '( concatenate ' &1,kstnm ' ' &1,kzdate ' ' &1,kztime ' ) '
title '( concatenate ' %ds ' ', r = '  ' %rs km ' ', az = ' ' %a2 ' ) '
*
xlabel "Time @(s@)"
ylabel "Displacement @(m@)"
*
ylim all
*
r $1$.z $1$.r $1$.t
*
taper w .2
bp bu c 0.02 0.067 n 4 p 1
qdp 1200
p1
qdp 1200
ylim off
title off
Contents


Example 2: Removing Instrument Responses

Here's another macro that manipulates some seismograms. This one reads in some seismograms, removes the instrument response uing a frequency-domain division (the transfer command), rotates the horizontal components into radial and transverse, and then cuts the signal using a minimum group velocity of 2 km/s. The instrument poles and zeros files are stored in the directory ../Resp (and are set up to produce displacements in meters). The original data are in the directory ../Usnsn and have the suffixes .LH[Z,E,N]. The macro takes two arguments: the first is a base file name for the three components; if the second argument is equal to all three components are processed. The arguments are referenced using $1, and $2 in the macro.










* * start with the vertical component * r ../Usnsn/$1$.LHZ * * set filter parameters for instrument deconvolution * setbb f1 0.005 setbb f2 0.006 setbb f3 45 setbb f4 55 * rmean taper transfer from polezero subtype ../Resp/$1$.LHZ.RESP freq %f1 %f2 %f3 %f4 w $1$_dsp.z * if $2 eq "all" * r ../Usnsn/$1$.LHN rmean taper transfer from polezero subtype ../Resp/$1$.LHN.RESP freq %f1 %f2 %f3 %f4 w $1$_dsp.n * r ../Usnsn/$1$.LHE rmean taper transfer from polezero subtype ../Resp/$1$.LHE.RESP freq %f1 %f2 %f3 %f4 w $1$_dsp.e * r $1$_dsp.n ch cmpaz 0 ch cmpinc 90 wh * r $1$_dsp.e ch cmpaz 90 ch cmpinc 90 wh * r $1$_dsp.n $1$_dsp.e rotate ch kcmpnm Radial w $1$_dsp.r $1$_dsp.t r $1$_dsp.t ch kcmpnm Trnsvrs wh * endif * * window from the origin time to OT + distance/2 * evaluate to t0 &1,o evaluate to t1 &1,dist / 2.0 + %t0 * cut %t0 %t1 r $1$_dsp.z $1$_dsp.r $1$_dsp.t * * convert from meters to cm * mul 100 * w over * cut off * * run a "system command" to remove the intermediate n & e files. * sc /bin/rm $1$_dsp.n $1$_dsp.e  
Contents



 






Example 3: Preparing Data for and Empirical Green's Function Analysis

Here's yet another macro that windows the data using a distance-dependent group velocity window, and sets up some of the origin time information in the SAC header. Aaron Velasco and I wrote this macro to prepare data for empirical Green's functions analyses. We designed the macro to be edited for each event that we were going to analyze. We store the group velocity window in the header, and used theose timess later to window the surface waves out of the original files. We also had a file naming convention that is easily changed.




 
$KEYS infile
*
*  macro to rotate and set headers
*   in the long period data
*
*  usage: macro lp_rotate infile base_sacfile_name
*
*  EDIT the Origin time variables: nyr, njday, ...
*       and the eqname variable.
*
*  if
*    files are name 9206281157.bks.?.lp
*  then
*     base_sacfile_name = 9206281157.bks
*
*  output:
*   sac files: KSTNM_lp_z KSTNM_lp_r KSTNM_lp_z
*   plot of the 3 components
*   group velocity windows are marked in the sac headers
*   for rayleigh (radial & vertical) and love waves (transverse)
*   to cut the surface wave window out use "cut t0 t1" on
*   the resulting files
*
* set up path to Original data
*
setbb orig /r2/thorne/gopher/nicaragua.9.2
*
*  set up the origin time
*
setbb nyr 1992
setbb njday 246
setbb nhr 00
setbb nmin 16
setbb nsec 00
setbb nmilsec 500
*
* Set up eq name and location
*
setbb eqname Nica_MS
setbb evlat 11.700
setbb evlon -87.300
*
* near and far distance (used for group velocity windows)
*
setbb n_f_dist 35.
*
* Rayleigh near and far velocity windows
*
setbb vrmin_f 3.2
setbb vrmax_f 4.4
setbb vrmin_n 2.15
setbb vrmax_n 4.4
*
* Love near and far velocity windows
*
setbb vgmin_f 3.8
setbb vgmax_f 5.3
setbb vgmin_n 2.4
setbb vgmax_n 5.3
*
*
qdp off
setbb ZC  %orig%/$infile$.z.lp
setbb NC  %orig%/$infile$.n.lp
setbb EC  %orig%/$infile$.e.lp
r %NC
setbb end1 (ADD &1,e -1.)
setbb beg1 (ADD &1,b +1.)
r %EC
setbb end2 (ADD &1,e -1.)
setbb beg2 (ADD &1,b +1.)
setbb emin (INT (MIN %end1 %end2%))
setbb bmax (INT (MAX %beg1 %beg2%))
cut %bmax%. %emin%.
r %NC %EC
cut off
ch evlo %evlon
ch evla %evlat
ch cmpinc 90
wh
setbb sta  &1,kstnm
rmean
rot
w %sta%_lp_r %sta%_lp_t
r %sta%_lp_r
ch kcmpnm radial
ch O gmt %nyr %njday %nhr %nmin %nsec %nmilsec
setbb evdist &1,gcarc
*
if %evdist gt %n_f_dist
markt v %vrmax_f %vrmin_f
endif
*
if %evdist lt %n_f_dist
markt v %vrmax_n %vrmin_n
endif
*
wh
r %sta%_lp_t
ch kcmpnm trnsvrse
ch O gmt %nyr %njday %nhr %nmin %nsec %nmilsec
*
if &1,gcarc gt %n_f_dist then
markt v %vgmax_f %vgmin_f
endif
*
if &1,gcarc lt %n_f_dist then
markt v %vgmax_n %vgmin_n
endif
*
wh
r %ZC
ch evlo %evlon
ch evla %evlat
rmean
ch kcmpnm vertical
ch O gmt %nyr %njday %nhr %nmin %nsec %nmilsec
*
if %evdist gt %n_f_dist
markt v %vrmax_f %vrmin_f
endif
*
if %evdist lt %n_f_dist
markt v %vrmax_n %vrmin_n
endif
*
w %sta%_lp_z
r %sta%_lp_z %sta%_lp_r %sta%_lp_t
ch kevnm %eqname
fileid l ul t list kstnm kcmpnm kzdate kztime
wh
synch
setbb shift &1,O
ch allt ( MUL -1 %shift )
title '( CONC 'Raw Velocity, Distance ' &1,gcarc ', Azimuth ' &1,az )'
p1
w over






Contents

Top | CJA's Guide Index | Ammon's Home | Department of Geosciences

Prepared by: Chuck Ammon
cammon@geosc.psu.edu
August 1996