Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
  • Apr 04 21:43
    jorcubillos commented #2594
  • Apr 04 20:37
    jorcubillos commented #2594
  • Apr 04 20:01
    jkmacc-LANL commented #2594
  • Apr 04 19:53
    jkmacc-LANL labeled #2594
  • Apr 04 19:53
    jkmacc-LANL assigned #2594
  • Apr 04 19:53
    jkmacc-LANL edited #2594
  • Apr 04 19:23
    jorcubillos opened #2594
  • Apr 03 13:51
    heavelock commented #2591
  • Apr 03 13:49
    WayneCrawford labeled #2593
  • Apr 03 13:48
    WayneCrawford assigned #2593
  • Apr 03 13:47
    WayneCrawford labeled #2593
  • Apr 03 12:09
    WayneCrawford commented #2593
  • Apr 03 10:36
    WayneCrawford review_requested #2593
  • Apr 03 10:35
    WayneCrawford opened #2593
  • Apr 03 10:20
    megies commented #2592
  • Apr 03 08:59
    megies commented #2591
  • Apr 02 17:03
    claudiodsf commented #2577
  • Apr 02 16:54
    aringler-usgs commented #2592
  • Apr 02 16:50
    amkearns-usgs synchronize #2592
  • Apr 02 16:00
    heavelock commented #2591
Stephen Hicks
@shicks-seismo
Screenshot 2020-03-25 at 14.19.31.png
Almost there ... but not quite. The colorbar is aligning with both "inner" subplots and then being overlain by the spectrogram
Thomas Lecocq
@ThomasLecocq
might be easier to call the specgram method of matplotlib youself directly :-)
ax.colorbar ?
and not fig ?
you can specify either which axes to steal space from or set up the axes for the colorbar yourself and use that
the "stealing space" approach might be hard to get right with a complex subplot layout, dunno. you could set up gridded subplots with 7/8 width for one column of axes and 1/8 of the width for tthe column of axes for the colorbars
dsentinel
@dsentinel
@ThomasLecocq, thanks. I'll get them started with that.
Stephen Hicks
@shicks-seismo
Thanks @megies & @ThomasLecocq - will try playing around with some of those ideas.
By the way, as a separate issues, I'm trying to vary the clip parameter in spectrogram and I see no apparent changes to how the spectrogram looks for vastly different lower and upper clip levels, both for dbscale=False and dbscale=True 😕😕😕
Stephen Hicks
@shicks-seismo
Screenshot 2020-03-26 at 14.18.50.png

Hurray, the subplot GridSpec configuration worked.
Here's the code snippet:

    fig = plt.figure(figsize=(12,14))
    outer = gridspec.GridSpec(len(st), 1, hspace=0.3)
    # Loop over trace from each channel in stream
    for n, tr in enumerate(st):
        inner = gridspec.GridSpecFromSubplotSpec(ncols=16, nrows=2, subplot_spec=outer[n], wspace=0.1, hspace=0.0)

        # Plot trace
        ax_trace = plt.Subplot(fig, inner[0,0:-1])
        ax_trace.semilogy(st_tr[n].times(), smooth(envelope(st_tr[n].data), 250), color='darkred')
        ax_trace.set_xlim([0, 3600])
        ax_trace.text(0.01, 0.95,'{:}.{:}'.format(tr.stats.station, tr.stats.channel),
                                            horizontalalignment='left', verticalalignment='top',
                                            transform = ax_trace.transAxes, fontsize=14)
        ax_trace.set_xticklabels([])
        ax_trace.set_ylabel("Amplitude")
        ax_trace.grid()
        fig.add_subplot(ax_trace)

        # Plot spectrogram
        ax_spec = plt.Subplot(fig, inner[1,0:-1])
        spectrogram(data=st_spec[n].data, samp_rate=st_spec[n].stats.sampling_rate, axes=ax_spec, show=False,
            dbscale=False, cmap="plasma", wlen=100, clip=[0.0, 1.0])
        ax_spec.set_ylabel("Frequency  (Hz)")
        ax_spec.set_xlabel("Time (s)")
        ax_spec.grid()
        fig.add_subplot(ax_spec)

        # Plot colorbar to spectrogram
        ax_cbar = plt.Subplot(fig, inner[1,-1])
        mappable = ax_spec.images[0]
        plt.colorbar(mappable=mappable, cax=ax_cbar)
        ax_cbar.set_ylabel("Spectral power")
        fig.add_subplot(ax_cbar)

We could always add something like this as an ObsPy function if anyone thinks it might be useful?

Tobias Megies
@megies
@shicks-seismo if you feel like it you could add it as a Tutorial code snippet http://docs.obspy.org/master/tutorial/index.html?highlight=tutorial
79seismo
@79seismo
Hi All, I am trying to extract some infor from my StationXML which contains information from two instruments (location codes 00 and 01) with each containing 3 channels (BH1, BH2, or BHN, BHE, or HHN, HHE, and BHZ or HHZ). First, I don't know which channel code is in the StationXML file, and because I am processing data from a large number of stations, I can't do the checking manually. I want to remove the instrument response by running tr.select for which I need to know the exact channel codes in the stationXML file. So, how do I extract these from the StationXML? Secondly, I want to read station/channel specific parameters: lat, lon, depth, dip, and azimuth from the StationXML and write them into SAC header. I tried to do this by xmlf = read_inventory(myXML); xmlf[0]. azimuth but that did not work. Many thanks!
Calum Chamberlain
@calum-chamberlain
Have you looked at the Inventory hierarchy? It you look at the figure and subsequent code blocks you should see how you access those attributes.
79seismo
@79seismo

Thanks @calum-chamberlain. Something like this xmlf[0][0][0] is what I am looking for, but still, the information is in a block (see below). Is it possible to separate different fields in the block? I tried a couple of things but it didn't work.

Channel 'BH1', Location '00'
Time range: 2018-07-06T20:00:00.000000Z - --
Latitude: 2.04, Longitude: -157.45, Elevation: 19.0 m, Local Depth: 1.0 m
Azimuth: 316.00 degrees from north, clockwise
Dip: 0.00 degrees down from horizontal
Channel types: CONTINUOUS, GEOPHYSICAL
Sampling Rate: 40.00 Hz
Sensor (Description): None (Trillium 360)
Response information available

Calum Chamberlain
@calum-chamberlain
Everything is an attribute of the objects. The Inventory is a list of Networkobjects, which in turn contains Station objects, which in turn contain Channel objects. Look at the Channel documentation.
What you are printing is the string representation of the Channel. Each of the things printed is it's own attribute of the Channel. For example, to get the channel-code, you would access the code attribute: a la: xmlf[0][0][0].code would be "BH1" in your example
79seismo
@79seismo
@calum-chamberlain Thank you! That is my problem. I don't exactly know the name of the attribute, e.g., if I use xmlf[0][0][0].Azimuth (which seems alright when I look at the output of xmlf[0][0][0]), it won't return anything. How do I figure out the exact name given to an attribute?
Calum Chamberlain
@calum-chamberlain

It's worth noting that all the Inventory, Network, Station and Channel objects are more than holders for lists and have other properties and methods as well - if you haven't tried already, hitting tab in iPython gives you a list of possible attributes.

As a general rule, only classes should have capital letters in them in Python (e.g. using PascalCase), whereas attributes, functions and methods should not have capitals and use underscores (e.g. snake_case). Obspy almost always follows that rule, so channel.Azimuth is actually channel.azimuth.

Nevetheless, the way to figure out the names is to read the docs. The list of attributes on the Channel docs tells you what they are called.
79seismo
@79seismo
Lovely! Many thanks @calum-chamberlain. Works like a charm!
Calum Chamberlain
@calum-chamberlain
Also, did you try the st.remove_response(inventory=xmlf) method - that should figure out the correct response for the correct data.
79seismo
@79seismo
yes, but I am not sure how I can extract the 6 channels separately to write SAC files after st.remove_response()
Calum Chamberlain
@calum-chamberlain
Ah okay, and you have more than just those 6 channels in your stream I guess, and you want some of the header information in the sac file. :+1:
79seismo
@79seismo
I only have either 3 or 6 channels in a single mseed from one station, but I need those attributes. Thanks a lot @calum-chamberlain .
Tobias Megies
@megies
@79seismo the Inventory class has methods for introspection / getting specific bits of information. Indexing your way through a big inventory is not what you want to do in 99% of cases. See http://docs.obspy.org/packages/autogen/obspy.core.inventory.inventory.Inventory.html, specifically select, get_orientation etc..
79seismo
@79seismo
@megies cool! Thanks!
Tobias Megies
@megies
In general, I'd either use these methods if you have some loops over your data and want to pick some piece of info out of a big inventory in that loop, or other use cases are looping over the inventory object (for net in inv: for sta in net: ...) and doing stuff automated in that loop. Depends on the workflow/application. But indexing should only happen after a fine grained select() to be 100% sure you pick the right piece of info. Even then you'd have to make sure you don't have multiple matching channels, to be super safe.
79seismo
@79seismo
I downloaded data tagged "best quality" from IRIS-DMC and so I am generally confident about mseed attributes. Thanks @megies !
79seismo
@79seismo
Hi, I am converting my traces into SAC files using st.write('fname', format='SAC') but having read docs, it is not clear to me how I should add headers like evla to my SAC file. Could you please explain? Thanks!
Thomas Lecocq
@ThomasLecocq
which docs did you read ?
Tobias Megies
@megies
@79seismo try reading a SAC file (e.g. st = read('/path/to/test.sac'))and look at st[0].stats.sac. You can set values there, but be aware you are entering the danger zone, you need to know yourself what should not be touched, especially the timing related variables in that header are really confusing and complex, you could easily mess up the timing info of your data
all of that is Trace based obviously
Unless you are bound to a specific existing workflow, you should really keep your waveforms separate from station/event metadata, though. There is a reason why all data centers keep it separate
http://docs.obspy.org/master/packages/autogen/obspy.io.sac.sactrace.html has some more info but that is a lower level interface, you're pretty much on your own there
79seismo
@79seismo
Thanks @megies does that mean I have to write the sac file, re-read it back on to Obspy and then write headers? @ThomasLecocq here's an example https://lists.swapbytes.de/archives/obspy-users/2016-February/001991.html
jkmacc-LANL
@jkmacc-LANL
If the Trace came from miniSEED, you are free to add a SAC header as a dictionary:
tr.stats.sac = dict(evla=32, evlo=105)
tr.write(..., format=‘SAC’)
79seismo
@79seismo
Thanks @jkmacc-LANL. My trace is the output of tr.remove_response(), and when I try tr.stats.sac = dict(evla=32, evlo=105), it returns the following error: AttributeError: 'Stream' object has no attribute 'stats'
jkmacc-LANL
@jkmacc-LANL
Your object is a Stream, not a Trace. Yes, you’ll need to add the SAC header to each Trace in a loop.
79seismo
@79seismo
@jkmacc-LANL I am a little puzzled because when I do len(tr) it returns 1, which means there's only one trace? How do I convert a stream into a trace? e.g. tr1 = Trace(tr) ? Any help? Thanks!
jkmacc-LANL
@jkmacc-LANL
A Stream is like a list of Traces. You would retrieve a Trace out of it just like you would retrieve an item from a list, even if the list contains only 1 item. In your case: tr = st[0].
79seismo
@79seismo
Yay! @jkmacc-LANL Thanks a lot!
jkmacc-LANL
@jkmacc-LANL
@79seismo This type of error is basic Python and ObsPy. The first line of the Stream docstring is List like object of multiple ObsPy Trace objects. I would encourage you to look into the documentation and the errors you’re getting before posting questions here. It will ultimately help you get your work done faster.
79seismo
@79seismo
@jkmacc-LANL appreciate it. I am slowly transitioning into Obspy from SAC (I am forced to due to the data format change at IRIS-DMC), and I usually post questions after attempting to solve problems a couple of times. In this instance, I did try out a few solutions. Python isn't the most intuitive language to me unlike fortran at least for now. Thanks for bearing with me.
jkmacc-LANL
@jkmacc-LANL
Please forgive my shortness; I hope you find the information you need to help you with your transition here and in the ObsPy docs. It is a welcoming community, and it has brought in researchers of many skill levels. This is a good thing, but it brings with it many questions, which can be a challenge. Good luck!
Try to name your variables accordingly (st / stream / st_... vs. tr / trace etc), it's for your own good and important especially if you share code with others.

Python isn't the most intuitive language to me unlike fortran at least for now.

I hope you don't want your GOTOs back. :stuck_out_tongue_winking_eye:

79seismo
@79seismo
Thanks @megies lets see!