CJA's Guide Index | Ammon's Home | Department of Geosciences
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
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
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
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 scalingContents
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
To set up a blackboard or temporary variable:
setbb my_variable 25
To access a header value of the first file
For the third file:
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)
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
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 offContents
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
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 overContents