Workflow examples - geoscience-community-codes/GISMO GitHub Wiki

Workflow example 1 In this example we:

  • load multiple event files (in this case, SAC files) into MATLAB as waveform objects
  • extract from a waveform object to create another waveform object
  • cross-correlate waveform objects
  • plot our aligned waveform objects

Workflow example 2 Load events from a Seisan database into a Catalog, and add waveform objects

  1. Create a variable that gives the path to the Seisan REA directory and the database where the S-files are stored: dbpath = '/raid/data/Montserrat/MontserratSeismicMastering/seisan/REA/MVOE_';

  2. Specify start and end dates/times (here we just load one day): startt = datenum(1996,11,1); endt = datenum(1996,11,2);

  3. Load the S-files into a Catalog object: cobj = Catalog.retrieve('seisan','dbpath',dbpath,'startTime', startt,'endTime', endt);

Note this actually creates a Seisan_Catalog object, which can do everything a Catalog object does, but has a few extra methods too.

  1. See the available properties - these are the variables that are part of a Seisan_Catalog object: properties(cobj)

You will see that this includes all the usual stuff like otime (origin time), lat, lon, depth, mag, mag_type (e.g. mb), etype (event type) & numberOfEvents. In addition, ontime and offtime store the beginning and end time of the waveform file belonging to that event. sfilepath stores the path to the S-file, wavfilepath stores the path(s) to the corresponding waveform files. arrivals store phase arrivals loaded from the S-files as GISMO Arrival objects. Since there will be many arrivals per event, and there are many events per Catalog, arrivals is a cell array (one element per event) of Arrival objects (each Arrival object can store all arrivals for one event). wavfilepath is similarly a cell array (1 element per event) of cell arrays (1 element per waveform file).

  1. See the available methods - these are built-in functions you can apply: methods(cobj)

There are about 30 methods in total. Here we will look at eev and addwaveforms.

  1. eev works in a similar way to the Seisan program of the same name. To start it, type:

cobj.eev()

or

eev(cobj)

These are functionally equivalent, but the first is in object-oriented style, the latter in the normal MATLAB style.

eev allows you to move through the list of events one at a time. Typing 'h' will show all options. 'S' will show the corresponding S-file. 'P' will plot the corresponding waveform file(s). You can jump to any time with t followed by the year, month, day, hour etc.

  1. addwaveforms will attach the waveform files to the Catalog. Again, there are two ways to do this:

cobj.addwaveforms()

or

addwaveforms(cobj)

Note that if you have a lot of events and/or a lot of channels, adding waveforms can be really slow.

  1. Now you can save the entire Catalog - including arrivals & waveforms, to MAT-files for fast loading next time:

save mycatalog.mat cobj

Note, this might be too slow or you might not have enough memory if your Catalog & set of waveform files is large, so break it up into pieces - the subset method is useful here, e.g. help Catalog/subset.

======

Additional, for the Seisan database from the Montserrat Volcano Observatory (MVO) only:

I worked for many years as the Seismic Network Manager at MVO and we added extra features into the MVO Seisan database. For each event we ran a program we wrote called ampengfft. This C-program computed amplitude, energy and frequency information for each channel in the Seisan waveform file, and added these measurements into the S-file as additional lines. This is the reason for the additional aef property in a Seisan_Catalog. I'm adding methods that use these data. So far these are:

regress_energy_vs_amplitude - determine relationship between max amplitude and total energy of a waveform for one event regress_energy_vs_amplitude_all - same but for a entire Catalog object

The reason for this investigation is that not all event types will scale in the same way as tectonic earthquakes. It depends on the envelope shape of the signal. Alas, it looks like I have the BGS version of the MVO database and not many files had ampengfft applied.