ObsPy: a bridge for seismology into the scientific Python ecosystem

Seismology is the study of the propagation of elastic waves through mostly solid media. The seismic waves result from natural (tectonic earthquakes, volcanoes, ocean waves) or man-made (nuclear explosions, quarry blasts, induced earthquakes) sources and are measured and recorded at seismographs. These perform point measurements of the elastic wavefield in up to three orthogonal directions. Each direction or component measures either the displacement, velocity or acceleration of the ground motion in the form of a one-dimensional, time-dependent signal. The resulting equally sampled time series are called seismograms or seismic waveforms.

With the emergence of digital seismology in the last couple of decades, numerous seismic waveform file formats have surfaced. Some, like the MiniSEED format [18] for data archiving and streaming, were designed with a specific purpose in mind, while the majority of formats are by-products of seismic signal processing or analysis packages that used custom formats for their I/O. The adoption of some of these software packages caused their I/O formats to become widely used data exchange formats—a purpose for which they were not designed. A recurring problem, for example, is the byte ordering: most formats do not specify whether they are written in big or little endian and for many formats, examples of both can be found. This myriad of file formats and software suites only accepting one particular file format caused many format converters to be written and distributed.

Each available waveform format in essence stores one or more time series and some meta-information about it. The essential ingredients of these formats can be distilled to a common subset of information enabling the use of a unified internal representation. ObsPy provides an implementation of that idea, which is described in the following section.

Assuming that a seismic data recording starts at time A, ends at time B and is equidistantly sampled without gaps and overlaps, it can be uniquely described by the time of the first sample, the sampling rate, the instrument it was recorded with and an array containing the signal. This is the common core information and all waveform file formats have to contain it in some fashion. In seismology receivers are uniquely identified by the so-called SEED identifiers [18] consisting of four short strings. A globally coordinated effort attempts to ensure that this system stays intact and as a direct result, most waveform formats adhere to it, offering a solution to the receiver location issue.

Seismic data also often have gaps and overlaps meaning that the data is not evenly sampled and the previous assumption does not hold true anymore. Gaps can result from short losses of power at the recording station or problems with the data transfer. Clock drifts and corrections can be one of the reason for the overlap of data points. Some data formats and ObsPy deal with this by allowing the storage of multiple chunks of waveform data, each being equally sampled and generally well behaved in itself.

Any additional information a waveform file format might support within ObsPy is stored separately as detailed in the next section.

Within ObsPy, waveform data are represented by a Stream object that acts as a container for any number of Trace objects. ObsPy defines a Trace to contain a single, contiguous, equally sampled time window of waveform data alongside the necessary meta-information. Each Trace object has a data attribute, which is a one-dimensional NumPy array. All further information is located in the dictionary-like stats attribute.

The stats object stores the four SEED identifiers network, station, location and channel, denoting the recording's physical location and instrument. It furthermore contains the times of the first and last sample, the data's sampling rate and interval and the number of samples. This information is partially redundant and ObsPy takes care that it stays consistent. For example, the endtime attribute is read-only and will be adjusted automatically if the start time, the sampling rate or the number of samples in the array has changed. Any additional information a specific file format might contain not covered by this abstraction will be stored in a container inside the stats object. See figure 1 for an illustration of the described internal representation.

Zoom In Zoom Out Reset image size

In order to support as many waveform file formats as feasible, a modular plug-in approach has been chosen with each supported data format being its own submodule and registering with ObsPy with the help of pkg-util's plug-in system. Listing 1 shows the registration of an example waveform plug-in. A file format submodule has to implement two functions with specified interfaces, and a third one is optional:

• is_format(filename): format identification function. Returns True if the passed file is of the submodule's format; False otherwise.
• read_format(filename, ∗∗kwargs): reads the file and returns a Stream object containing a representation of the file's data.
• write_format(stream, filename, ∗∗kwargs): writes a Stream object to the given filename; this function is optional, and if not given, no write support for the format will be available.

The forced modularity of this strategy fosters a clean separation of concerns with each submodule being implemented and tested independently. In addition, it allows users to extend ObsPy with I/O support for formats that are not used widely enough to justify integration into the main ObsPy library. A common example for this is output formats of numerical waveform solvers. pkg-utils' plug-in system will register these additional formats for a seamless integration with the rest of ObsPy.

ObsPy comes with a top-level read() function, a single entry point when reading waveform data. See the listing 2 for a usage example.

Listing 2:.  snapshot of an interactive Python session demonstrating the usage of ObsPy's read() function. It will detect the file's format and call the appropriate reading routine. In case the read() routine detects a valid HTTP URL it will download the resource before proceeding. In case it detects an archive format it will be decompressed first.

>>> import obspy

>>> st=obspy.read("filename")

>>> st

$\lt$ obspy.core.stream.Stream at 0x2f284d0 $\gt$

>>> print st

3 Trace(s) in Stream:

BW.RJOB..EHZ $|$ 2009--08-24T00:20:03.000000 Z -... $|$ 100.0 Hz, 3000 samples

BW.RJOB..EHN $|$ 2009--08-24T00:20:03.000000 Z -... $|$ 100.0 Hz, 3000 samples

BW.RJOB..EHE $|$ 2009--08-24T00:20:03.000000 Z -... $|$ 100.0 Hz, 3000 samples

The read(filename, ∗∗kwargs) routine calls all registered formats' is_fileformat() functions until one returns True. Depending on the format in question the is_fileformat() routine parses the first couple of bytes or performs more complicated heuristics. For performance reasons the format detection routines have to be fast. On top of that, ObsPy performs the format detection according to a manually curated list so that the most commonly used ones are tested first, improving average performance. After the format has been determined the appropriate format's read_fileformat() function will be called, which parses the file and returns a Stream object. The format detection can also be skipped by providing the format to the read() routine.

This structure permits the sharing of capabilities among all formats by implementing them as part of the read() routine. Examples of this are automatic file downloading if a valid HTTP URL is detected and the decompression of a number of different archive formats.

Web services provided by data centers around the globe are an important source for waveform data. ObsPy implements clients able to interact with a comprehensive selection of such data centers [20]. A waveform data request will end up being stored in a Stream object so the workflow following the data acquisition is identical no matter the origin of the data.

A unified internal data representation opens the possibility of defining methods transforming it. With respect to this, ObsPy offers a comprehensive collection of signal processing routines frequently used in seismology by relying heavily on functionality coming with NumPy and SciPy. These are implemented to match the needs seismologists have for a certain processing operation while keeping the data self-consistent. An example of this is the Trace.decimate() method which will apply a filter, decimate the data and adjust the sampling rate metadata of the times series. Caused by the potentially large number of samples, most operations are implemented as in-place modifications of the Stream and Trace objects.

Most processing methods are defined on a single Trace; the same methods on the Stream objects will call the corresponding method on each of their Trace children. This allows for a natural handling of three and more component data. A number of methods are specific to Stream methods. Consult table 1 for a selection of available signal processing routines.

In order to offer a large variety of functionality, ObsPy once again employs a plug-in system. Listings 3–5 demonstrate this by example with the Trace.taper() method. The employed plug-in system empowers the usage of functionality defined by different modules into a simple domain-specific API usable by scientists. It furthermore encourages code deduplication and reuse and thus enables ObsPy to offer a large variety of different functions without having to implement and test the details in many cases.

The calculation of seismic travel times is a problem frequently encountered in seismology. If an earthquake occurs at point X the question is what seismic phases arrive at what time at point Y. The cost for full 3D simulations of the whole Earth [29] are excessive for many applications and less accurate but faster approximations are often sufficient. In case a 1D spherically symmetric Earth model is applicable, the frequently used method of Buland and Chapman [3] enables very fast travel time calculations. It has been implemented in a Fortran package commonly called iaspei-tau [14, 26]. A more feature-rich Java version of this algorithm is also available [4]. The obspy.taup submodule uses ctypes from the Python standard library to provide a pythonic interface to iaspei-tau; see the listing 6. Figure 2 plots the travel times of an earthquake through the Earth versus the geographic distance in degrees.

Zoom In Zoom Out Reset image size

Seismic receivers distributed around the globe aim to measure and record ground motion as accurately as possible. Many factors influence the characteristics of the final waveform, among them the frequency response of the physical instrument, the effects of any amplifiers, of analog and digital filters and of the digitalization. For many applications studying the Earth, it is crucial to eliminate these effects to get the best possible estimate of the true ground motion. The first step when correcting data for the influence of the seismic receiver and processing chain is to calculate the frequency response of the recording system. The seismograms are then deconvolved with this response to obtain a seismogram with physical units unbiased by instrumental effects in the frequency band of consideration. This process is known as instrument correction in seismology.

Seismic recording systems are described by a linear chain of different stages or elements. The SEED data format [18] accurately represents these, is widely accepted by the community and has been in use since the early 1990s. Recently, StationXML [19], the designated successor of SEED, was developed, and the community is starting to adopt it. Being an XML format, it has many benefits in comparison to the binary SEED format regarding tool support, ease of data distribution and human readability. Except for some minor differences, both formats store the same information.

The standard workflow for calculating the frequency response of a seismic receiver system from SEED files involves rdseed [12] to generate RESP files, a textual strict subset of SEED. These files are then fed into evalresp [11] resulting in the frequency response. The detour via the RESP files is necessary as it is the only input file format evalresp accepts. Performing an instrument correction with metadata stored in StationXML files involves one additional step—the conversion of the StationXML files to SEED files using yet another tool.

Conducting this in Python involves a number of system calls and unnecessary I/O operations making it ill-suited to modern large data workflows. Replacing evalresp would be a major effort as many pitfalls exist when calculating instrument responses, which can greatly influence the final result. Instead ObsPy offers a direct bridge from StationXML (parsed through lxml) to evalresp utilizing ctypes to call internal evalresp functions. This approach is not feasible for code that is still actively developed as no guarantees on the stability of the internal API are usually granted. However, evalresp is not in active development anymore and only receives an occasional maintenance commit.

ObsPy's internal representation of StationXML data is used to derive the nested C structures evalresp's internal functions expect. These are defined and initialized using ctypes as demonstrated in listings 7 and 8. Ctypes also requires a definition of the function headers as shown in listings 9 and 10. Figure 3 displays the amplitude and phase response calculated from a StationXML file also illustrating the effect of different stages in the recording system.

Zoom In Zoom Out Reset image size

A seismic receiver's recording system usually has several stages resulting in a fairly complex representation within evalresp. Evalresp's file parsing stage takes care of translating the SEED constructs to the appropriate internal structure and sometimes makes non-obvious decisions. To ensure the correctness of the evalresp integration we performed an extensive test [16]. We define our integration of evalresp acting on StationXML files to be correct if the final response is equivalent to a response calculated by converting StationXML to SEED and SEED to RESP files on which evalresp is acting. We downloaded almost the complete set of StationXML inventory data available from IRIS, the largest data distributor in seismology. Data from over 27 000 stations, most with multiple recording channels defined for different periods in time result in well over 100 000 instrument responses. We reiterated our implementation until the tests passed, giving confidence that our solution can be safely applied to any data encountered in the world wide community.

Contrary to the well-established, virtually bug-free legacy codes described above, many scientists using Python in their analytic workflow tend to develop bits and pieces of codes that are applied to their input data, and thus at first sight not useful for others. This general ascertainment is also true in seismology where input data can have different formats, outputs need to be compatible with different downstream processing software and so forth. In this section, we first present a case study of the translation of a C code that can read an exotic format called Win and then present a research code developed during a PhD thesis that is now integrated in ObsPy as a new module. Thanks to the pluggability of ObsPy, only the necessary processing steps need to be written.

WIN is the format imagined by a Japanese company named Hakusan as default storage for their Datamark datalogger series. The data within a one-minute WIN file is highly compressed, every second being compressed with a different level from the previous. Up to now, there were three ways to convert WIN to a more standard format, namely SAC: Win2Sac and japan2sac in Linux and Win2Sac GUI on Windows. Win2Sac codes come from Japan, while the japan2sac was written as part of a research project in Chambery (France) and is not completely finished. In order to process these data using state of the art routines and software, e.g. using ObsPy-based scripts, one has to convert the whole archive to SAC, which is readable by ObsPy. Each step of this process doubles the volume of the archives stored on disk, and potentially adds errors to the final product.

Thankfully, sources of Win2Sac are available online [22] and together with the datasheet [30] describing the WIN format, we have been able to write an ObsPy plug-in to read it. The final goal was to be able to read directly from the archive without duplication and process those data the same way as data originating from other file formats. Even if one's goal is to convert WIN to MiniSEED, the best solution today will be to write a Python script using ObsPy and no longer rely on a succession of steps.

The state of a volcano and its unrests can be evidenced by the amount of amplitude (or energy) recorded by seismic sensors on its flanks. Real-time seismic amplitude measurement (RSAM (1)) was initially introduced by Endo and Murray [8] to forecast eruptions and assess the state of volcanic activity. One can also compute the standard deviation of the squared amplitudes, or real-time seismic energy measurement (RSEM (2)) [5], of the signals over a certain time interval divided by the number of measurements.

$$\begin{eqnarray}&&{\rm RSAM}=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{n}(|{{A}_{i}}|),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm RSAM}=\sqrt{\frac{1}{N}\mathop{\sum }\limits_{i=1}^{n}{{({{A}_{i}}-{{A}_{{\rm avg}}})}^{2}}},\end{eqnarray} \tag{ 2 }$$

where A_i is the seismic signal's amplitude, A_avg is the average over N samples.

As multiple sources are overlapping, volcano-seismologists typically calculate the RSEM or RSAM values for separate, previously defined frequency bands (spectral seismic energy measurement (SSEM) [28], similar to SSAM introduced by Stephens [27]. The data are first demeaned and bandpass filtered with a Butterworth filter of order two. Each value is then calculated on a 30 s window (100 × 30 samples). Finally, the daily 10th and 25th percentiles and the median are typically computed to partially remove undesired transients and earthquakes which are considered a disturbance of the tremor data [7].

The RSAM, RSEM, SSAM and SSEM algorithms have been implemented in the obspy.signal.ssxm submodule. The final code takes advantage of the resampling abilities of Pandas (Python Data Analysis Library), resulting in a total of less then 50 lines of code for the ssxm() method.

4.3.1. Example results

The interest of filtering the signal in different frequency bands can be illustrated with the Kawah Ijen volcano (East Java, Indonesia) seismic data. The seismic time series are generally polluted by anthropic activities during the day, especially in high frequency bands ($\gt$10 Hz) and particularly for the stations located closer to the source of 'cultural' noise, in this case the parking and base camp for climbing the volcano. One may even observe weakened working activity during the prayer day (i.e. Friday, see figure 4). To minimize the effects of transients, one can use the 10th percentile (p10) of the amplitude in a specific frequency band (figure 5). Except during volcanic unrest (figure 5), the p10 data from station PSG are still highly affected by this anthropic noise, compared to POS and PUN, which are located around the crater, thus further away from the human activity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Other studies investigated some particular processes by filtering the seismic signal. For example, SSEM computations might provide a measure of the rate of strain released, if persistent fracturing at any scale produces a sustained seismic signal [6]. One can also assess the stationarity of the microseism noise sources (0.1–1 Hz) to evaluate the results from the velocity variations using ambient seismic noise cross correlation techniques [17].