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()
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
next_fast_len
when computing ffts
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
).
times
is a function, so st[0].times()
is the way to get the values
get_events
method to catalog (e.g. use event = catalog.get_events(...)
) that supports the same arguments as the FDSN client.get_events
method,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
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?
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.
Stream
object are exactly the same
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
.
:(
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.
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:
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