Hi everyone,
I have been experimenting with the PPSD spectral_estimation.py code, adding some colorization to the various plot overlays such as the mean, mode, and percentile graphs. I want to add a little bit of optimization to the PPSD plot, as the overlays default to black. How would I go about becoming a first-time contributor to the project ?
Something weird happens when I try to trim a trace (both with trim and slide commands). The start-time of the trimmed output is like 4 milliseconds ahead of the start-time that I provide in the code, irrespective of the time I give (unless I give an integer value in second+ say 4 millisecond). I kinda need to trim different stations' data all at the same point in time to run a cross-correlation code on it... but this glitch is keeping me from it.
I tried playing around with resampling, but that doesn't work.
Anybody else face such issues?? How did you solve it??
nearest_sample that sample might be left or right of the given time. If you need to have the same time and the samples dont align you could resample your data to have the exact same sampling points, but personally I would not do that and rather adjust the timing in the cross correlation function that you get as outputhttps://docs.obspy.org/master/packages/autogen/obspy.core.trace.Trace.trim.html?highlight=nearest_sample
I am an obspy beginner. I want to get the events catalog from ISC, including arrival information and pick-trace information, and then use this information to download the corresponding waveform from IRIS. (I am labeling the dataset for deep learning training).
I am encountering some problems here,
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
-->client = Client("ISC")
-->t0 = UTCDateTime(2020, 1, 2)
-->t1 = UTCDateTime(2020, 1, 4)
-->cat = client.get_events(t0, t1, minlatitude=-30, maxlatitude=-16,
minlongitude=20, maxlongitude=28, includearrivals=True)
-->print(cat.str(print_all=True))
-->print(cat.events[0].picks[3])
The output I am getting has an error in the Network code. WaveformStreamID
network_code: 'IR'
station_code: 'KHWEE'
channel_code: ''
location_code: '--'
The network is supposed to be BX here. Now I cannot get waveforms from IRIS, I will be happy to discuss this further because I am getting lots of errors. Maybe because I am a newbie.
from obspy.clients.fdsn import Client
client = Client('IRIS')
inv = client.get_stations(network='BX', level='station')
sta_codes_bx = [sta.code for net in inv for sta in net]
for pick in event.picks:
is_IR = pick.waveform_id.network_code == 'IR'
should_be_BX = pick.waveform_id.station_code in sta_codes_bx
if is_IR and should_be_BX:
pick.waveform_id.network_code = 'BX'
client.get_waveforms(....)
Hi everyone,
If I use the command scan inside a script.py it raises an error like this:
PermissionError: [WinError 5] Access denied: '/$Recycle.Bin\S-1-5-18'
on Ubuntu the error is almost the same.
What I do is the following:
from obspy.imaging.scripts.scan import scan
scan("path_to_traces","output_path",verbose=True)
Trying to emulate what I do from the shell:
obspy-scan input_trace_path -v -o outputTracePath
From the shell I do not have any problem, from the script it doesn't work.
Any clue?
Thanks
6 replies
I want to make a recod_section. According to the page obspy.core.stream.Stream.plot in the "section" part, I must define distance or event coordinates to plot section type. Is there a documentation or example code you would know that can help me how do I add them in the sac header?
3 replies
6 replies
Hello all,
I need help with instrument response removal using poles and zeros. All the information I have about the seismometer is stated here (https://www.gfz-potsdam.de/en/section/geophysical-imaging/infrastructure/geophysical-instrument-pool-potsdam-gipp/pool-components/poles-and-zeros/3d-geophone) and below:
No Damping Resistor
Gs = 28.8 V/(m/s); Ds = 0.56; fs = 4.5 Hz
2 Zeros
z1 = 0
z2 = 0
2 Poles
p1 = -15.88 23.43i
p2 = -15.88 -23.43i
Normalisation Factor
a = 1.0
I can see the code example using PAZ dictionary here: https://docs.obspy.org/tutorial/code_snippets/seismometer_correction_simulation.html#using-a-paz-dictionary but not sure which of these parameters correspond to the ones used in the example.
There is also a response file available here: https://www.gfz-potsdam.de/fileadmin/gfz/sec22/pdf_doc/GIPP/RESP.XX.NS532..SLZ.HG6.375.100000.4_5B.28_8.txt and also on IRIS NRL: http://ds.iris.edu/NRL/sensors/hgs/RESP.XX.NS532..SLZ.HG6.375.9090.4_5B.27_7 but didn't work with ObsPy.
I kept getting "ValueError: No matching response information found" error when I use these files (after modifying the channel name, station name and start and end dates).
Kindly provide me some guidance on using the parameters or response file successfully.
Thanks!
6 replies
4 replies
if key == 'NOTE':
try:
value = [cleanup_and_decode_string(line)
for line in value.split(self.line_terminator)
if line]
except:
value = ''
else:
value = cleanup_and_decode_string(value)
Can I save the network information in individual file and write it to all files?
3 replies
Hello, I'm trying to process seismic data from MASW using the SWPROCESS library, it intrinsically uses some obspy libraries. Following the tutorial on the page, which loads a series of .dat data in SEG2 format separated by offset, I am not successful in inverting my data. I've already tried to use my data in SEG2, SEGY and SU format and still no success. With that, my last attempt was to try with the SEG2 format named the same as the one in the tutorial. However, I get the following error:
# Path (relative or full) to a folder containing the data files. Data files must be in either the SEG2 and/or SU data format.
path_to_folder = "data/wghs/"
set0 = [f"{path_to_folder}{x}.dat" for x in range(6, 7)]
fnames_set = [set0]
print(fnames_set)
for set_number, fnames in enumerate(fnames_set):
print(f" set{set_number} includes {len(fnames)} files from {fnames[0]} to {fnames[-1]}")
names = None
workflow = "time-domain"
transform = "fdbf"
fmin, fmax = 3, 100
vmin, vmax, nvel, vspace = 100, 500, 400, "linear"
fdbf_weighting = "sqrt"
fdbf_steering = "cylindrical"
snr = True
noise_begin, noise_end = -0.5, 0.
signal_begin, signal_end = 0., 0.5
%matplotlib qt5
settings = swprocess.Masw.create_settings_dict(workflow=workflow,
trim=trim, trim_begin=trim_begin, trim_end=trim_end,
mute=mute, method=method, window_kwargs=window_kwargs,
transform=transform, fmin=fmin, fmax=fmax, pad=pad, df=df,
vmin=vmin, vmax=vmax, nvel=nvel, vspace=vspace,
weighting=fdbf_weighting, steering=fdbf_steering,
snr=snr, noise_begin=noise_begin, noise_end=noise_end,
signal_begin=signal_begin, signal_end=signal_end,
pad_snr = pad_snr, df_snr=df_snr)
start = time.perf_counter()
wavefieldtransforms = []
for fnames in fnames_set:
wavefieldtransforms.append(swprocess.Masw.run(fnames=fnames, settings=settings))
print(wavefieldtransforms)
end = time.perf_counter()
print(f"Elapsed Time (s): {round(end-start,2)}")From that I get the following error:
Could anyone help me to solve this problem? I would be extremely grateful!
2 replies
I am trying to analyse daily seismic recordings near a river to capture signals associated with river processes. It is new for me, and not sure if I am doing things right. I cut the signal in 6 hours windows and tried to plot spectrogram/PSD (not sure of the difference between Spectro and PSD).
Here is the basic python script I used.
import numpy as np
import matplotlib.pyplot as plt
import obspy
from obspy import read
st = read("WU.JSP2.00.ELZ.D.2022.231.000002.SAC_6h")
print(st)
tr=st[0]
print(tr.stats)
sps =int(st[0].stats.sampling_rate)
st.plot()
fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2]) #[left bottom width height]
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])
t = np.arange(tr.stats.npts) / tr.stats.sampling_rate
ax1.plot(t, tr.copy().data, 'k')
tr.spectrogram(wlen =2sps, per_lap=0.95, dbscale=True, log=True,cmap="rainbow")
st = read("WU.JSP2.00.ELZ.D.2022.231.000002.SAC_12h")
print(st)
tr=st[0]
print(tr.stats)
sps =int(st[0].stats.sampling_rate)
st.plot()
fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2]) #[left bottom width height]
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])
t = np.arange(tr.stats.npts) / tr.stats.sampling_rate
ax1.plot(t, tr.copy().data, 'k')
tr.spectrogram(wlen =2sps, per_lap=0.95, dbscale=True, log=True,cmap="rainbow")
I don't know how to include the colour scale too. Any advice will be great.
Thx.
We would like to inform you about a public beta release of a new FDSNWS station service at the Italian EIDA node (@INGV). This new services implementation should remove some important performance limitation of the currently used version. The service is available at the following URL: https://webservices.ingv.it/beta/fdsnws/station/1/.
We also would very much appreciate any feedback, bug reports or other comments, in order to catch any issue before the production release occurs which is planned within the next few weeks. Thanks!
1 reply
