--------SAC software as illustrated by the 2001 Anza earthquake-------- On October 31, 2001 (Halloween!), a M=5.1 earthquake occurred near Anza which was widely felt across San Diego County. Here is the event info from the Southern California Seismic Network (SCSN) operated by the USGS and Caltech: 9718013 le 2001/10/31,07:56:16.630 33.5083 -116.5143 15.22 5.09 l 1.0 The first number is the "cusp id" which is a unique number that is assigned to each event located by the SCSN. Confusingly, these numbers are not in any kind of recognizable order. The date and origin time (UT) are given by 2001/10/31,07:56:16.630 and then followed by the location and depth in km. The 5.09 gives the assigned magnitude. I'm not sure what the "le" and the "l 1.0" are. (local magnitude?). There is also information about this quake on a website (http://eqinfo.ucsd.edu/special_events/2001/304/a/index.html) maintained by the Anza group here at IGPP. Waveform data for this event may be found on the Macs in ~shearer/PROG/CLASS/227C/SAC/ANZA_5.1/9718013 and on the Suns in /net/chip/227C/ANZA_5.1/9718013 The contents of this directory look like: 9718013.AZ.BZN.BHE.sac 9718013.CI.FUR.BHN.sac 9718013.CI.RRX.BHN.sac 9718013.AZ.BZN.BHN.sac 9718013.CI.FUR.BHZ.sac 9718013.CI.RRX.BHZ.sac 9718013.AZ.BZN.BHZ.sac 9718013.CI.FUR.HHE.sac 9718013.CI.RRX.HHE.sac 9718013.AZ.BZN.HHE.sac 9718013.CI.FUR.HHN.sac 9718013.CI.RRX.HHN.sac 9718013.AZ.BZN.HHN.sac 9718013.CI.FUR.HHZ.sac 9718013.CI.RRX.HHZ.sac 9718013.AZ.BZN.HHZ.sac 9718013.CI.FUR.HLE.sac 9718013.CI.RRX.HLE.sac 9718013.AZ.CRY.BHE.sac 9718013.CI.FUR.HLN.sac 9718013.CI.RRX.HLN.sac 9718013.AZ.CRY.BHN.sac 9718013.CI.FUR.HLZ.sac 9718013.CI.RRX.HLZ.sac etc. Each of these files is a single time series from a seismograph that recorded this event. These data may be readily obtained from the SCEC data center (www.scecdc.scec.org) by using the STP program (http://www.data.scec.org/STP/stp.html). The "standalone" version of the program is most convenient--it downloads the waveforms directly to our local machines. On the Suns, my copy of the program is located at ~shearer/SCEC/STP/stp. On the Macs, you can download STP from Kris Walker's website at: http://sail.ucsd.edu/~walker/sac.tar.gz You will then have to move the sac directory into /usr/local on your machine ("sudo mv sac /usr/local" may be necessary to have permission to do this) and add setenv SACAUX /usr/local/sac/aux to your .cshrc file. NOTE: If you are running SAC on grunt in the Barnyard, then you will need to use: setenv SACAUX /usr/local/lib/SAC/aux because the directory structure is a little different. The data are in a standard format called "SAC" as indicated by the .sac suffix on each file name (more about SAC later). The file naming convention is as follows: 9718013 = cuspid AZ = network (in this case, IGPP's very own Anza network) BZN = station name B = channel B = broadband H = highly broadband (seems basically the same as B but with a higher sampling rate) V = short period E = also short period (V and E seem interchangable) A = accelerometer H = gain (H = highgain, L = lowgain, S = accelerometer) E = component (E=east, N=north, Z=vertical) The SAC format was created by seismologists at Lawrence Livermore National Laboratory (LLNL) (my hometown!) in the 1980s as part of a seismic analysis package. There is a SAC manual for the 1988 version in the bookcase in the Barnyard. Documentation for the latest version is available on the web at www.llnl.gov/sac/ and the program is also now being distributed by IRIS at http://www.iris.edu/manuals/sac/ You can also obtain help within the program by entering "help command" where command is the name of the command that you would like explained. SAC has become the most widely used format for storing seismic data. Alternatives include the "ah" format developed at Lamont and GFS (Guy's filing system) used by many here at IGPP. We have conversion programs that convert SAC and AH files to the GFS format that may be helpful to you in seeing how the formats work if you are interested. The SAC program is installed on the Suns. To run it, you will need to modify your .cshrc file to include: setenv SACAUX /fish/local/src/sac/aux setenv SACDIR /fish/local/src/sac/bin EXAMPLES OF SAC COMMANDS I am a novice at SAC, but have managed to get a few things to work. Here are some examples. First, go to the directory where the data are, in this case ~shearer/PROG/CLASS/227C/SAC/ANZA_5.1/9718013, then enter sac: shearer@katmai 56> sac SEISMIC ANALYSIS CODE [03/01/2005 (Version 100.00)] Copyright 1995 Regents of the University of California SAC> To plot data, first open a plot window with the bg (begin graphics) command: SAC> bg x This will open an x-windows window on your screen for the plot output. ---------------------------------------------------------------- NOTE: You may get the following error message: SAC> bg x > ERROR 203: Can't create X window. Check DISPLAY environmental. In this case, enter: setenv DISPLAY :0.0 before running SAC and you should be OK. ----------------------------------------------------------------- Next, read in a seismogram: SAC> read 9718013.AZ.BZN.BHZ.sac Next, plot the seismogram: SAC> plot (NOTE: On the Suns, this does not work if Netscape is also running!) This station is very close to the event and is severely clipped for this large event. We can list the header information in the SAC file using the "list header" command: SAC> listhdr FILE: 9718013.AZ.BZN.BHZ.sac - 1 ---------------------------- NPTS = 10218 B = -2.123000e+01 E = 2.341950e+02 IFTYPE = TIME SERIES FILE LEVEN = TRUE DELTA = 2.500000e-02 IDEP = UNKNOWN DEPMIN = -8.388608e+06 DEPMAX = 8.202399e+06 DEPMEN = 1.691010e+05 OMARKER = 0 KZDATE = OCT 31 (304), 2001 KZTIME = 07:56:16.630 KINST = BHZ KSTNM = BZN CMPAZ = 0.000000e+00 CMPINC = 0.000000e+00 STLA = 3.349150e+01 STLO = -1.166670e+02 STEL = 1.301000e+03 STDP = 0.000000e+00 EVLA = 3.350834e+01 EVLO = -1.165143e+02 EVDP = 1.522000e+01 DIST = 1.430736e+01 AZ = 2.625446e+02 BAZ = 8.246018e+01 GCARC = 1.286557e-01 LOVROK = TRUE NVHDR = 6 LPSPOL = FALSE LCALDA = TRUE KCMPNM = BHZ KNETWK = AZ Note that the source-receiver range (DIST) is only 14 km. Anza station CRY is not much further away (22 km), but its waveform is much less severely clipped, e.g. SAC> read 9718013.AZ.CRY.BHZ.sac SAC> plot Note that each time "read" is called in these examples, it replaces any files already in memory. To get around this, use the "more" option, e.g. SAC> read 9718013.AZ.BZN.BHZ.sac SAC> read more 9718013.AZ.CRY.BHZ.sac SAC> plot1 In this case the "plot1" option is used to plot both traces at once. The default time window for the plots will be the length of the time series. To change this, use the xlim command: SAC> xlim 0 50 SAC> plot The plot will now be between 0 and 50 s. You can also truncate data when they are read into memory using the "cut" command: SAC> cut -10 50 SAC> read 9718013.AZ.CRY.BHZ.sac SAC> plot If xlim has not been set, the plot will default to between -10 and 50. Data outside this window have not been read into memory and cannot be plotted even by changing xlim. The SAC program defaults to what it calls "quick and dirty" plots that skip points so that they go faster on the screen. This can be deceiving because vital details (e.g., high frequencies!) in the plot can be lost. Thus, I recommend that you turn off this feature by entering: SAC> qdp off before plotting. On modern workstations, the plots are still plenty fast. At 100 km, the station on Mount Soledad provides a nice P and S wave: SAC> read 9718013.AZ.SOL.BHZ.sac SAC> plot This time series can be filtered by using the bandpass command: SAC> bandpass corner 1 5 SAC> plot This filters between 1 and 5 Hz. Subsequent filters will operate on this filtered time series, not the original. To restore the original, just enter read again with no argument. There are also "highpass" and "lowpass" and "demean" commands. For example, SAC> read SAC> highpass corner 5 SAC> plot will plot the Soledad record above 5 Hz. It would get tedious if every file name had to be read in separately. Fortunately, SAC allows the use of wildcards to read in multiple file names all at once. For example, SAC> read *BHZ* will read in all of the BHZ data (Warning: this will be lots of traces). The plot command will then plot them one at a time--you hit 'return' to move to the next trace. If you don't want to plot the remaining traces, then enter 'k' or 'kill'. If you want to plot them all without stopping, enter 'g' or 'go' (but this goes so fast on my computer as to be pretty useless). To read in all of the ANZA BHZ data: SAC> read *AZ*BHZ* Note that multiple wildcards can be used. ? can be used as a wildcard for a single character. For more about wildcard options, enter "help wild" Filtering and other processing commands will operate on all of the time series in memory at once. Jeanne finds this a convenient way to filter data. The filtered SAC files may be written out as, for example, SAC> write temp.sac Making a hardcopy of your plot is rather awkward. Jeanne tells me the following will work: SAC> begindevice sgf (your window will disappear) SAC> plot1 (or whatever...) This will create a file called f001.sgf (and then f002.sgf, etc if you make further plots). To translate into postscript: % sgftops f001.sgf f001.ps This will of course have problems if you don't have write permission into the directory that you are running SAC from (often the data directory). To get around this, use the "read dir" option to read the data from a different directory (see below). To plot multiple traces vs. range on the same plot (this is called a seismic "record section"), you must enter the Signal Stacking Subprocess: SAC> read *AZ*BHZ* SAC> sss Signal Stacking Subprocess. SAC/SSS> plotrecordsection A reduction velocity of 8 km/s can be invoked as follows: SAC/SSS> vm 1 refractedwave vapp 8 SAC/SSS> plotrecordsection The time window can be changed as follows: SAC/SSS> timewindow 20 50 SAC/SSS> plotrecordsection To leave the subprocess, use the QUITSUB command: SAC/SSS> quitsub SAC> To quit SAC, use the QUIT command: SAC> quit EXAMPLE: PLOTTING RECORD SECTION OF ANZA EARTHQUAKE I could not find a nice map of the Tri-Net stations with station names on the map. However, from the Tri-Net web site I was able to find maps that would identify the stations when I clicked on them. From this, I made up a list of stations extending north from the earthquake. To read these data in, it is most convenient to use a macro file so that changes can be easily made. This file is called "sac.macro1" (in ~shearer/PROG/CLASS/227C/SAC/) and looks like this: read dir /net/chip/ANZ/9718013.full *MSJ.BHZ* read more *SVD.BHZ* read more *SBPX.BHZ* read more *LUG.BHZ* read more *VTV.BHZ* read more *LKL.BHZ* read more *EDW.BHZ* read more *LRL.BHZ* read more *JRC.BHZ* read more *CWC.BHZ* read more *TIN.BHZ* Note the use of the "read dir" option to change the directory to read the data from to the specified directory rather than the current directory (the one that you started the SAC program from). This remains in effect for the subsequent reads and is a convenient option so that one can save macros files, etc. in the current directory. NOTE: The Anza stations have a different sample rate than the SCSN/TriNet stations. This causes problems with the 'sss' command in SAC (see below). For this reason, I include only SCSN/TriNet stations in the above macro. (This problem could be avoided if the traces were resampled using the 'interpolate' command. For example, to resample to 100 Hz, enter 'interpolate delta 0.01'. I have not yet tried this to see how well it would work in our example....) Notice also the use of the more command to avoid overwriting the previous data. These stations are in order of increasing range. Within SAC, one can simple enter: SAC> macro sac.macro1 to read the data, provided this file is in your current directory. Otherwise, enter the complete path name or use the SETMACRO command to set a default directory for the macro commands. Next, the traces can be plotted all at once using plot1: SAC> plot1 If desired, they can be filtered: SAC> lowpass corner 1 SAC> plot1 to see frequencies below 1 Hz. A true record section plot can be produced, for example, as follows: SAC> sss Signal Stacking Subprocess. SAC/SSS> plotrecordsection SAC/SSS> timewindow 0 200 SAC/SSS> plotrecordsection SAC/SSS> vm 1 refractedwave vapp 8 SAC/SSS> plotrecordsection The 'ls' command, when used within SAC/SSS, lists the station ranges which can be convenient for station selection. OTHER SAC COMMANDS rotate This can be used to rotate horizontal component data into radial and transverse components. I have been having trouble getting this to work.... help Lists all the SAC commmands. Then do "help command" to learn more about each one. ylim Change y plotting limits