.. _sac1: ********************************* Seismic Analysis Code 1 ********************************* One tool that geophysicists use to process seismic data is Seismic Analysis Code (SAC). SAC is installed on all of the computers in the Mac Lab, and can be obtained for no charge from IRIS (Incorporated Research Institutions for Seismology, which is one of the main consortiums for archiving and distributing seismic data) if you are affiliated with an IRIS institution. SAC has many seismic processing functions, makes several kinds of plots, and can read and write ".sac" data files (more on this later). SAC has been around for a while, which means it has a few limitations based on its age, but it is still used enough that geophysicists should be familiar with it. SAC capabilities include many tools commonly used when processing seismic time series data: arithmetic, Fourier transforms, integration/differentiation, spectral estimation, filtering, stacking, decimation and interpolation, correlation, phase picking, instrument correction, vector component rotation, envelope creation, frequency/wave number analysis, and plotting. While we will not be covering all of these in this class, if you have a working knowledge of how to use SAC, you can figure out how to use any of these tools without too much trouble. There are many alternatives that replicate the functionality of SAC, so while you should know how to use SAC, you may want to explore some of these other options if your research involves processing seismic data. First, there is a MATLAB toolbox MatSeis that has many functions useful for seismic data analysis. There are also are many individual MATLAB functions that people have written to perform many of the operations in SAC. You can of course also write your own MATLAB programs to process seismic data. As a Python fan, I like using the Python package ObsPy to do many of the same things that SAC does (I believe that ObsPy is installed on the computers in the Mac Lab). And for serious data processing, there are many C and Fortran codes available to perform various tasks. Talk to other graduate students to find out how they process data. ================ Running SAC ================ Unlike MATLAB, SAC does not have a nice GUI. You run SAC from the Unix shell by typing ``sac``, followed by a carriage return. The Terminal should now look something like this: :: SEISMIC ANALYSIS CODE [11/11/2013 (Version 101.6a)] Copyright 1995 Regents of the University of California SAC> We now have a command prompt within SAC to enter commands. To quit SAC, type either ``exit`` or ``quit``, followed by a carriage return. This will return you to the Terminal prompt. We will be reading and writing files from within SAC. By default, SAC reads and writes files from the current working directory. You will probably want to create a separate folder for all of your SAC work in this class. One method is to create a directory ``SAC`` in the ``Documents`` folder, though you are welcome to arrange things however you like. ================= SAC Basics ================= SAC is used to analyze time series data. Like MATLAB, SAC is interactive and command-driven (either by entering commands directly in the terminal, or by executing commands saved in a script), so you need to give it commands for it to do anything. SAC commands can be typed either as upper case or lower case (commands are all converted to upper case before they are interpreted by SAC). However, if you are giving SAC file names to read from or write to, SAC interprets the case of those statements literally. There are also abbreviated names for many functions to save you time in typing them. I will generally use all uppercase and give the full function name when introducing a function, but lowercase when making subsequent references to it. I will mention the abbreviated names as well, but will stick to using the full names in the notes. SAC commands are fall into three different categories: (1) parameter-setting, which changes the values of the parameters within SAC, (2) action-producing, which carry out actions based on these parameter values, and (3) data-related, which determine the data in memory that these actions will be performed on. Note that this means that commands that set parameters appear to do nothing until you execute the corresponding action-producing command. Additionally, there is the ``HELP`` command, which lets us read the documentation on any of these commands. ``help`` without any additional arguments brings up a list of all commands. ``help `` brings up the documentation page for that particular command. To exit the documentation and return to the SAC command line, type ``q``. ---------- FUNCGEN ---------- To put a basic time series into SAC, use the ``FUNCGEN`` command. This can be entered in all caps, in all lowercase, or with the shorthand ``fg``. If you simply type ``funcgen``, followed by a carriage return, SAC will create a default time series for you. There are additional options that can be specified (we will get to those in a second), but this is a basic way to create data within SAC. After you type ``funcgen``, it is not clear that anything has happened. To see the results of the ``funcgen`` command, we need to be able to plot the data and look at it. --------- PLOT --------- To plot the current data that is in memory, we use the ``PLOT`` command, or ``p`` for short. Try this with the function that you created using ``funcgen`` -- the plot should show up in an external window, and look something like this: .. image:: sac1_impulse.png To change the limits on a plot, use the ``XLIM`` and ``YLIM`` commands. The syntax is ``xlim ``, then issue the plot command (remember, SAC commands to set parameters do nothing until you enter a corresponding action-producing command). Also note that this will affect all future calls to ``plot``, even if you generate new data in SAC. To turn off axis limits and return to the default behavior, enter the command ``xlim off``. --------------- Some Details --------------- Now that we can create functions and plot data to visualize, here is a bit more about what is going on. When we invoke ``funcgen`` with the default parameters, SAC creates an impulse function (also known as a delta function) in memory. By default, time ranges from 0 to 99 with 100 data points, and the impulse occurs at :math:`{t=50}`. The function is one at the time of the impulse, and zero otherwise. When we plot the data, SAC creates a plot in an external window. The plots are relatively simple, and reflect the legacy of SAC being originally used on Tektronix 401X graphics terminals: .. image:: sac1_tektronix.png These 1980's era terminals were state-of-the-art when SAC was developed at Lawrence Livermore National Laboratory, and so some limitations of the graphics remain in the present version of SAC -- in particular, SAC will often desample a time series to plot it more quickly, which was a necessity on the Tektronix system, but not a real concern on modern computers. However, we will see that many of these limitations have seen improvements in the modern version of SAC. ----------------- FUNCGEN (again) ----------------- Now that we know how to view our results, we can play a bit more with the ``funcgen`` command. Look at ``help funcgen`` to see some of the additional options. In addition to the default ``impulse`` function, we can also generate ``step``, ``boxcar``, ``triangle``, ``sine``, ``line``, ``quadratic``, ``cubic``, ``seismogram``, ``random``, and ``impstrin`` (multiple impulses). Some of these have additional parameters; for example, for ``sine`` you specify two parameters (the frequency and phase), for example: :: SAC> funcgen sine 0.05 0 SAC> plot will produce the following graph: .. image:: sac1_sine.png Other parameters that you can specify when calling ``funcgen`` (regardless of the type of function you are creating) include the number of points ``npts``, the time resolution ``delta``, and the beginning time ``begin``. The default values of these parameters are ``npts = 100``, ``delta = 1``, and ``begin = 0``. To call ``funcgen`` with additional parameters, use ``funcgen impulse npts 100 delta 1 begin 0``. There are also default values for all of the functions that take parameters. Once you set a parameter to a different value, that becomes the default value for subsequent calls to ``funcgen``, unless you specify that parameter again. **Exercise:** Try out some of the other functions available in SAC, in particular some of the functions that require additional parameters like ``sine`` and ``random``. Also produce functions over differing time window using ``npts``, ``delta``, and ``begin``. Make a plot of each function to verify that you know how ``funcgen`` works. ================================ READ, WRITE, and ".sac" Files ================================ One thing you may have noticed is that when you create a new function in memory, you lose all ability to access the previous function. However, modern seismology often requires processing multiple signals simultaneously. How can we do this in SAC? In order to get multiple signals into memory in SAC, we need to read data from a file using the ``READ`` command (``r`` in shorthand). SAC has its own data format, and is one of the many standard data formats that you will find for waveform data. Other formats include "ah" (Ad-Hoc, the data format of the IRIS/PASSCAL program), "SEED" (Standard for the Exchange of Earthquake Data, standard format of the IRIS Data Management Center), "SUDS" (Seismic Uniform Data System, developed by the USGS), "SEG-Y" (the standard for seismic reflection data), and others that continue to pop up from time to time. There are additional software packages for many of these additional formats that fulfill many of the same functions as SAC. SAC files all have a header that displays information about the data in the file. Once you read in the data, you can access the header information using the ``LISTHDR`` command, which can be abbreviated with ``lh``. The header for a seismic data file will usually include such information as time ranges, instrument details, and any seismic phase pick information that is present in the file. We can also look at the header of files that we create within SAC -- these typically only include time information and the type of function, but are still useful if you don't remember the details of how you created the data that you currently have loaded in memory. However, what if we don't have our data in a file, but want to generate it within SAC? Then we use the ``WRITE`` command (shortened to ``w``). Once we have created our function, we call ``write ``, and SAC writes the function to a .sac file with the specified name. Then we can read the file back into memory using ``read ``. When you call ``write`` with the filename of an existing file, it overwrites that file without a warning, so you need to be careful when using ``write``. **Exercise:** Generate a function using ``funcgen``, write it to disk using ``write``, and then read it back into SAC using ``read``, verifying that you have succeeded by plotting the data. Why go to the trouble of writing the signal to disk? Now that we can read and write data, we can get multiple signals into SAC using the ``read`` command, followed by multiple filenames: ``read file1.sac file2.sac file3.sac``. This will load all three signals into memory, allowing you to process and plot the signals. Alternatively, you can load files one at a time without overwriting the data presently in memory using ``read more `` (you can call ``read more`` with no files loaded in memory, so it is okay to call it without having already some files in memory). **Exercise:** Create three functions using ``funcgen``: an impulse, a boxcar, and a triangle. Save each one to disk, and then read all three into SAC. Plot the functions. Notice that when you use ``plot`` with multiple functions in memory, it plots them one at a time, waiting for you to hit return to see the next one. ==================================================== Obtaining ".sac" Files from Public Data Sources ==================================================== While you can learn to use SAC with data generated within the program, it is far more interesting to use real data from real earthquakes. One such source if data is the IRIS Data Management Center, which has a nice web interface called "Wilber 3" for finding earthquakes and downloading waveforms. Go to the Wilber 3 website at http://ds.iris.edu/wilber3/find_event to take a look at recent earthquake recordings that are available to download. Click on an event from the map or the list (by default, it shows events for the past 30 days, but there is a drop-down menu above the map that allows you to change what is displayed). From there, it will go to a map showing the epicenter and the location of stations. By default, it shows only broadband channels from the Global Seismic Network, but you can change the options at the right of the map. Below the map, you can select the stations that you want to download using the checkboxes at the left of the table. Click on "Request Data" to download the selected stations (you will need to select the desired time windows and data format; choose "SAC binary (little endian)" for use with SAC on the Mac Lab computers). You should receive an email link that allows you do download the data. **Exercise:** Once you have some real waveforms, try loading them into SAC and plotting them as above. You may find that it is a pain to type in the full filename when reading the data, but SAC recognizes the same wildcards as the shell, so you can use wildcards to reduce the amount of typing that you do when loading multiple files into memory. ================================ Plotting Multiple Functions ================================ What if we want to plot several functions on the same graph? SAC has additional plotting commands that handle multiple types of data. ``PLOT1`` (``p1`` when shortened) plots multiple signals on separate axes with a common time axis for all signals. If the three signals have different time ranges, SAC extends the time axis to be large enough to show all signals. Each graph has its own unique amplitude limits. However, if you would like to give each plot its own time axis, use the ``relative`` command (as opposed to the default ``absolute``). If you have three data files in memory and type ``plot1 relative`` into the command terminal, then it will give each plot a separate time axis that corresponds to the limits of the data. We can use ``PLOT2`` (``p2``) to plot multiple signals on the same axes, and the time and amplitude limits are chosen to be large enough to show all points. As with ``plot1``, you can have the plot show relative values instead of absolute values. With multiple plots (particularly ``plot2``), we may want to have our plots displayed in color. To do this, we need to use the ``color`` command. ``color`` has several options. If we just want to change the color of the lines in the plot, use the ``color `` command (see the ``help color`` documentation to the available colors). This will make all future plot calls use this color for the lines on the graph. If you want different colors for each successive line that you plot, then use ``color on increment on``, which cycles the colors for each line that is plotted. **Exercise:** Use the seismograms that you just downloaded to make plots using ``plot1`` and ``plot2``. First, read in all of the data files into SAC. Use the ``listhdr`` command to examine the header variables, which includes information on the instrument type, event, and the recorded signal. Plot the signals using both ``plot1`` and ``plot2``. You will note that it is difficult to distinguish the traces on the overlay plot, so change your ``plot2`` to be a color plot to help you distinguish the waveforms. Notice that in your data files, we have three components per station, which required several data files. This highlights one of the limitations of SAC -- each file can only contain one trace, so multiplexed data needs to be stored in multiple files. This is one of the downsides of the .sac format, a particular problem given the amount of data available today. Fortunately, automating things with a computer makes it easier to handle many files, so with appropriate care we can handle this efficiently and reliably and continue to take advantage of the many functions available in SAC.