ObsPy integration

[liamtoney]

ObsPy is the de facto standard tool for processing seismological data using Python. GMT is already heavily used by seismologists. It seems natural to make PyGMT integrate better with ObsPy. This issue can serve as a platform for broad discussion of how we might accomplish this.

On the seismology side in general, we already have Figure.meca() wrapped. Some GMT plotting commands (namely, sac) could be wrapped by PyGMT and integrated with standard ObsPy data objects, such as Trace and Stream objects. We could even configure Figure.plot() to take Inventory objects.

For waveform plotting, imagine:

import pygmt
import obspy

tr = obspy.read()[0]

fig = pygmt.Figure()
fig.coast(...)  # Make your map
fig.plot(...)  # Plot your stations / event / etc.
fig.sac(trace=tr, ...)  # Plot your waveform
fig.show()

(Of course, perhaps we'd want to rename Figure.sac() for this example...)

We'd want to make ObsPy an optional dependency if we went down this route. I think that most of the work would be related to I/O, type-checking, etc. — as well as surveying the seismology community on what'd be useful.

Do folks have other ideas? Pinging @megies and @krischer in case you have any thoughts on this (thanks).

[core-man]

Some GMT plotting commands (namely, sac) could be wrapped by PyGMT

As a seismologist, I almost use gmt sac in every project. Currently, I call lib.call_module to use gmt sac. See a simple example:

import pygmt
from pygmt.clib import Session

fig = pygmt.Figure()
fig.basemap(region=[-50, 100, -1, 3], 
            projection="X10c/5c",
            frame=['xa20f10+l"T (sec)"', 'ya1+l"NO."', "WSrt"])
with Session() as lib:
    lib.call_module("sac", "{} -En -Tt1 -C -M1.0c -W0.5p".format("saclist"))
fig.show()

saclist is a sac file list with the format:

filename [X Y [pen]]

$ cat saclist              
seis1.sac
seis2.sac
seis3.sac

If gmt sac can be wrapped, the script to plotting seismic waveforms will be simpler.

integrated with standard ObsPy data objects, such as Trace and Stream objects.

gmt sac can only plot sac format, but it seems that other formats, e.g., miniseed, become more and more and more popoular. We can still use gmt sac to plot those formats if something like I/O and type-checking can be well-resolved by pygmt via ObsPy, which supports a lot of seismic waveform format (maybe all). See Waveform Import/Export Plug-ins of ObsPy.

[liamtoney]

We can still use gmt sac to plot those formats if something like I/O and type-checking can be well-resolved by pygmt via ObsPy, which supports a lot of seismic waveform format (maybe all). See Waveform Import/Export Plug-ins of ObsPy.

This is a good point. If ObsPy is made a [optional] dependency, then we could leverage its I/O to read in all sorts of files (as well as Trace and Stream objects).

[core-man]

if gmt sac can be wrapped, the seismological community will benefit a lot from PyGMT and GMT.

@liamtoney Could this issue be completed in two or more PRs? For example,

1. Wrap gmt sac

   I think this is the first step before ObsPy can be a [optional] dependency. Meanwhile, seismologists could use
   PyGMT to plot SAC data a simple way before ObsPy is made as a dependency.

2. Plot seismic data in other formats via ObsPy I/O

   It may be a little difficult for PyGMT to support all the available waveform formats in ObsPy. So we have to think about
   which format PyGMT will have to support.

   Miniseed is one of the most popular formats in global seismology now, and miniseed (waveform data) + StationXML
   (metadata) is a good combination for a seismologist. ObsPy reads miniseed data to Trace/Stream objects, and StationXML metadata to be Inventory objects. Based on my experience, we should be able to use PyGMT to plot miniseed + StationXML if gmt sac is wrapped.

[liamtoney]

@core-man yes, exactly! The PyGMT devs were thinking that two PRs would make sense.

Agreed on your first one. For number 2, I'm not quite sure what you mean by:

It may be a little difficult for PyGMT to support all the available waveform formats in ObsPy.

I was assuming we'd just utilize the ObsPy I/O here. We could either read it to e.g. a NumPy array, or quietly change it to a SAC file and plot it that way (hacky, perhaps).

[seisman]

I was assuming we'd just utilize the ObsPy I/O here. We could either read it to e.g. a NumPy array, or quietly change it to a SAC file and plot it that way (hacky, perhaps).

The major difficulty is that other waveform formats (e.g., miniseed) don't have time markers like SAC does (e.g., T0 to T9 in SAC headers), so we can't use gmt sac -T option unless someone is clever to find the solution, but I agree we should support more formats.

[core-man]

I was assuming we'd just utilize the ObsPy I/O here. We could either read it to e.g. a NumPy array, or quietly change it to a SAC file and plot it that way (hacky, perhaps).

The major difficulty is that other waveform formats (e.g., miniseed) don't have time markers like SAC does (e.g., T0 to T9 in SAC headers), so we can't use gmt sac -T option unless someone is clever to find the solution, but I agree we should support more formats.

@liamtoney The above difficulty is also what I am concerned, while the format transformation itself should be easy as your suggestion. We don't have SAC headers related to time picks if we only use miniseed and StationXML, unless we also send an additional input with time picks to pygmt.

[willschlitzer]

I'm not familiar with SAC data/files, but is it information that could be relatively easily transferred to/from a numpy array/pandas dataframe? I figure sac has to be wrapped to accept this before we can have ObsPy integration?

[michaelgrund]

I'm not familiar with SAC data/files, but is it information that could be relatively easily transferred to/from a numpy array/pandas dataframe? I figure sac has to be wrapped to accept this before we can have ObsPy integration?

In principle everything from a SAC file could be transferred to a numpy array/pandas dataframe or any structure we need. However, some of these steps are already handled by ObsPy, thus using these resources potentialy could help us a lot. As Liam suggested, first we should discuss if:

• sac has to be wrapped as is with support of data only in SAC format or
• if we want to make ObsPy an optional/additional dependency and we simply use all available resources which would allow to support all data formats that ObsPy can work with