Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
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
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:
  1. 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,
  2. 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:

  1. Request download for time window a (for argument's sake let's say 1980-1990).
  2. The obspy Mass Downloader downloads waveforms and the station response with the channel information for the requested time window .
  3. 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