obspy/obspy

ObsPy: A Python Toolbox for seismology/seismological observatories.

Tobias Megies

@megies

i think it would really be good to have, but probably it would be good to have a discussion how to best do it or have some different options to do it, so I think we should try to get as many people working on this into the discussion (i know there are more people doing it)

igmuhlmann

@igmuhlmann

Jep, the issue will be the best platform. Didn't know about gist yet... ;)

geophysics91

@geophysics91

atleast can anybody suggest how to read ascii data, while reading above file i am getting error unknown file format

Thomas Lecocq

@ThomasLecocq

@geophysics91 https://earthscience.stackexchange.com/questions/18361/how-to-read-and-plot-seismic-waveform-data-in-ascii-format-using-obspy

geophysics91

@geophysics91

Thanks @ThomasLecocq .

Thomas Lecocq

@ThomasLecocq

make sure to create a full Trace object with all .stats etc.. and then set your data in .data

Pascal Audet

@paudetseis

Hi ObsPy developers - I have found myself having to write a custom method to shift a Trace object along the time axis for various applications (e.g., stacking along move-out curves, etc.). I'm wondering if this would be a good method to add to the Trace class? Here's the hack I've coded up to do this - let me know if this would deserve a proper PR (with tests):

import types
import numpy as np
from obspy import read

def shift(self, tshift):
    """
    Method to shift a :class:`~obspy.core.Trace` object along the time axis

    The `starttime` attribute is adjusted to preserve absolute timing. This
    method can be useful in applications where stacking along moveout curves
    is required. By using the FFT, this method correctly interpolates the
    amplitude between two time samples.

    :param tshift: time shift in seconds
    :type tshift: float

    """

    # Define frequencies
    faxis = np.fft.fftfreq(self.stats.npts, d=self.stats.delta)

    # Fourier transform
    ftrace = np.fft.fft(self.data)

    # Shift
    for i, freq in enumerate(faxis):
        ftrace[i] = ftrace[i]*np.exp(-2.*np.pi*1j*freq*tshift)

    # Back Fourier transform
    self.data = np.real(np.fft.ifft(ftrace))

    # Update stats.starttime
    self.stats.starttime -= tshift

if __name__ == "__main__":

    # Example usage of the `shift` method
    trace = read()[0]
    trace.plot()

    # Extend the trace instance with new method
    trace.shift = types.MethodType(shift, trace)

    # Apply a shift to a trace
    trace.shift(5.)
    trace.plot()

Thomas Lecocq

@ThomasLecocq

good idea, I wrote something similar for msnoise (for correcting timestamping of samples not aligned on delta grid) https://github.com/ROBelgium/MSNoise/blob/master/msnoise/api.py#L1770-L1804 -- Question: you're doing np.exp(-2.*np.pi...) and thus applying a negative shift for a positive tshift arg, (confirmed by self.stats.starttime -= tshift ... intuitively, I would shift "right" if tshift is positive

and you also might want to use next_fast_len when computing ffts

Pascal Audet

@paudetseis

Thanks for the feedback. Yes, the negative factor is just how I've been implementing things in my codes, so that a positive tshift moves the seismogram to the right. It might be more appropriate/intuitive to do the opposite, in which case the negative factors become positive (in the exp function and starttime).

Thomas Lecocq

@ThomasLecocq

intuitively , if I shift by -5 seconds, I expect my starttime to go back in time by 5 seconds,... but it's my Belgian interpretation :)

zeynepodabass

@zeynepodabass

Hi, when I take take data by using ObsPy.client I can't have a plot. <Figure size 800x250 with 1 Axes> just this occurs. Also when I try to see st[0].times it gives <bound method Trace.times of <obspy.core.trace.Trace object at 0x7fa1f9384820>> how can I fix these?

Thomas Lecocq

@ThomasLecocq

depending on how you call the plot function (and where) , you sometimes need to "show" it using matplotlib (plt.show)

and times is a function, so st[0].times() is the way to get the values

zeynepodabass

@zeynepodabass

I was using the same method but I couldn't have the figure. I am working on Jupyter maybe it is the problem. When I run the cell with %matplotlib command I have the figures.

79seismo

@79seismo

Hi, what does it mean when you define pre_filt=None in Trace.remove_response? Thanks!

Thomas Lecocq

@ThomasLecocq

That the raw data is NOT filtered before removing the response

spri902

@spri902

I'd like to pad my traces with 50 zeros at the start of each trace. I'm using cross-correlation for pick refinement in an automated phase picking algorithm and the issue I'm running into is that in some traces the p-arrival comes in early, at sample 5 or 10 for example, on other traces the p-arrival may come in at sample 100. I'm using a fixed window size in my cross-correlation routine (let's say 25 samples) and I'd like to add these zeros to the beginning of the signal to prevent an error because the window size is too large(i.e. the window tries to look backward 25 samples but there are only 5 or 10 because the arrival is so early). I realize I can use the trim method for this, but am unclear about usage. I've seen some comments about the method creating a masked array. What is that? Additionally, if I add these zeros to the front of each trace does it change the time vector?

79seismo

@79seismo

@ThomasLecocq Is that actually possible? My guess is some sort of filtering is needed to keep the deconvolution stable, e.g. spectral holes will distort the signal?

WING CHING

@jwjeremy

Hi i would like to ask are there any methods or functions that extract specific event from catalog object base on the event.resource_id? thanks !

Calum Chamberlain

@calum-chamberlain

@jwjeremy I don't think there is an Obspy native way to do this. Two simple options are to:

Install ObsPlus which adds a few features including a get_events method to catalog (e.g. use event = catalog.get_events(...)) that supports the same arguments as the FDSN client.get_events method,
If you just want a way to get events by event-id, the way I usually do this is to make a dictionary of my catalog a la:
```
cat_dict = {ev.resource_id.id: ev for ev in catalog}
```
then get events by key, e.g. event = cat_dict[event_id]

@spri902 you can just add the required zeros to the tr.data attribute. This will not change the tr.stats.starttime though, so you will need to reset that (e.g. tr.stats.starttime -= pad_length / tr.stats.sampling_rate where pad_length is the number of zeros pre-prended).

Bear in mind that adding zeros will affect your correlations, so if your preferred shift includes correlation of some of your waveform with zeros your correlation will be reduced.

RE masked array, see the numpy Masked Array docs

spri902

@spri902

@calum-chamberlain thanks! I'd like to pad with values that have the noise characteristics of my signal rather than zeros. (bc zeros running into my signal create an issue for my autopicker.) Any thoughts on pre-padding with "noise"?

Calum Chamberlain

@calum-chamberlain

@spri902 pre-padding with noise is definitely not a bad idea - I haven't done it personally, but you could try taking your data, computing the FFT, randomising the phase and doing the IFFT to get phase randomised noise with similar amplitude characteristics.

spri902

@spri902

@calum-chamberlain :) thanks

Thomas Lecocq

@ThomasLecocq

alternatively, you could just cut the first n seconds of your data, flip it left-right and append it at the beginning

vivigb

@vivigb

When we download the waveform data using obspy.clients.fdsn.mass_downloader, from IRIS , The response is already removed from it or we have to do it afterwards?

Thomas Lecocq

@ThomasLecocq

You have to do it afterwards

vivigb

@vivigb

So there is no optional method to mass-download the data with its response attached, I will get the data in the form of velocity/displacement etc. instead of counts?

Thomas Lecocq

@ThomasLecocq

you will get the data in counts/volts and you can tell mass_downloader to download the responses too (saved to stations/ by default I think), then you're free to do what you want to process the data , filtering, removing instrument response using your parameters etc

Stephen Hicks

@shicks-seismo

Yep, as @ThomasLecocq says, each station will have its own stationxml file saved in stations/ (the other directory being waveforms/ by default).
You can then use obspy.read_inventory followed by trace.remove_response to correct for the response of each station.

vivigb

@vivigb

Thank you so much! I will try that as you guys suggested.

vivigb

@vivigb

I downloaded some waveform files and corresponding response files from IRIS
Now I am trying to remove the response through a loop.Eventhough my attempt is successful.
I am getting an error for a file

Traceback (most recent call last):

  File "<ipython-input-72-9c07bc7f593d>", line 8, in <module>
    water_level=60)

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\stream.py", line 3155, in remove_response
    tr.remove_response(*args, **kwargs)

  File "<decorator-gen-147>", line 2, in remove_response

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\trace.py", line 273, in _add_processing_info
    result = func(*args, **kwargs)

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\trace.py", line 2885, in remove_response
    output=output, **kwargs)

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\inventory\response.py", line 1674, in get_evalresp_response
    freqs, output=output, start_stage=start_stage, end_stage=end_stage)

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\inventory\response.py", line 1632, in get_evalresp_response_for_frequencies
    end_stage=end_stage)

  File "C:\ProgramData\Anaconda3\envs\roses\lib\site-packages\obspy\core\inventory\response.py", line 1581, in _call_eval_resp_for_frequencies
    raise e('check_channel: ' + m)

ValueError: check_channel: Illegal RESP format


 EVRESP ERROR (... [File: <stdin>; Start date: ; Stage: 3]):
    check_channel; units mismatch between stages,
    skipping to next response now
 EVRESP ERROR (... [File: <stdin>; Start date: ; Stage: 3]):
    check_channel; units mismatch between stages,
    skipping to next response now
 EVRESP ERROR (... [File: <stdin>; Start date: ; Stage: 3]):
    check_channel; units mismatch between stages,
    skipping to next response now
 EVRESP ERROR (... [File: <stdin>; Start date: ; Stage: 3]):
    check_channel; units mismatch between stages,
    skipping to next response now
 EVRESP ERROR (... [File: <stdin>; Start date: ; Stage: 3]):
    check_channel; units mismatch between stages,
    skipping to next response now

Whats this mean?

16 replies

rameesmir

@rameesmir

Hello experts, I am using ppsd.calculate_histogram() to customize the data for PPSD estimation which works fine for ppsd.plot(). Apparently, the same isn't carried to ppsd.plot_spectrogram(). I use Obspy 1.1.1 with Python3.5 on Ubuntu 16.04 LTS. Any help will be deeply appreciated.

4 replies

vivigb

@vivigb

Hii everyone,
So even though obspy can read the SAC format traces, Obspy won't support removing the response from the sac format data , if the response file is in .PZ format?
or simply obspy won't support .PZ format response files?

Thomas Lecocq

@ThomasLecocq

I think all this works properly, did you try with a response file that WORKS ? (== not your Kazakstan buggy response xml as source ?)

ObsPy doesn't care about data formats, once it's read from SAC or from MSEED, the internals of the Stream object are exactly the same

vivigb

@vivigb

I think all this works properly, did you try with a response file that WORKS ? (== not your Kazakstan buggy response xml as source ?)

st=read('BNDL.BHZ.2016.151-0000.01SPS.SAC',attach_response='BNDL.BHZ.PZ')
st.remove_response(output = 'VEL')

I am trying to remove the response of sac file using .PZ response file.
But I think I have to give the inventory object. Since it is giving the following error :
ValueError: No response information found. Use "inventory" parameter to specify an inventory with response information.
But .PZ format is an unknown format for read_inventory.
:(

Emmanuel D. Castillo Taborda

@ecastillot

Good morning everyone!, I have been trying to initialize a FDSN Client with this link: http://168.176.34.103:18080/fdsnws but returns the next exception: obspy.clients.fdsn.header.FDSNException: No FDSN services could be discovered at 'http://168.176.34.103:18080/fdsnws'... I don't understand why it happens? can you help me please? thanks.

Thomas Lecocq

@ThomasLecocq

you can try to remove the trailing /fdsnws

or add a /

I've had issues with my SeisComP3 system, similatr to this

Emmanuel D. Castillo Taborda

@ecastillot

thanks a lot @ThomasLecocq ! It was literally like you said

Thomas Lecocq

@ThomasLecocq

(which one? )

2 replies

uky-jps

@uky-jps

Question, I am reading in several dataless seed files using obspy.io.xseed into individual Parser objects. Is there a way to combine them? I want to produce a single *.xml file for all the response information but I can't find any information on how to combine different parser objects into one before I write out the xml file. Thanks in advance.

Fred Massin

@FMassin

You can use new += other with station metadata inventories...

1 reply

Fred Massin

@FMassin

I think you need inventory objects, from read_inventory(), but I m not sure to understand your objective

1 reply

Fred Massin

@FMassin

Cool! There’s a quick tutorial about instrument metadata management, it might worth reading it ;)

Joel D. Simon

@joelsimon

Hi everyone -- first off, this is my first question and I've been using ObsPy
for a while so that's a testament to its cleanliness and documentation. Thanks
to all the contributors!

Here's my problem: I'm writing mseed files that have a time correction
applied. It is not obvious to me how to attach time correction info to the
metadata (e.g., like in trace.Stats()). I see that I can update the bool "Time
correction applied" activity header AFTER writing the mseed with
util.set_flags_in_fixed_headers(). However, my understanding is that this only
overwrites flags and does not also set their associated values (in my case,
bytes 40-43 in the mseed header hold the LONG of the actual time correction).
I'm playing with a short function to f.seek(40,0) the mseed file and then
struct.pack() the time correction value, but I want to ensure that I'm not
reinventing the wheel before I properly implement and test it.

Thank you.

1 reply

Peter Makus

@PeterMakus

Hi everyone,
First of all, thanks a lot for this amazing project. I've been using Obspy for about two years now and thanks to your really thorough documentation I haven't had any major issues yet.
I've got a particular question/request concerning the obspy.clients.fdsn.mass_downloader.

Here is what I would like to do:
I am creating a global database of waveforms and station response files (i.e., station XMLS) to be used in a receiver function imaging workflow. I want to be able to update my database at any time, but I encountered the following issue:

Request download for time window a (for argument's sake let's say 1980-1990).
The obspy Mass Downloader downloads waveforms and the station response with the channel information for the requested time window .
Request download for time window b (i.e., a different time window than a, e.g., 1995-2005).
1. The Obspy Mass downloader will again download the waveforms for time window and download a station response with the channel information for time window b. However, the old station XML will be deleted before! (This is actually the expected behaviour according to the Mass Downloader's documentation.)

Now, I would like to pass a function to stationxml_path that either updates the station xml (I want to keep the response information for both time window a or b!) or download the whole station xml for all available time windows every time I start a download for any time window (which would be the less elegant solution to my problem, but should be ok as well). What would the function have to return in order for this to work?

Thanks in advance,
Cheers,
Peter

6 replies

Where communities thrive

People

Repo info

Activity