Getting started question.
I think I imported obspy into my new anaconda install base environment
(base) PS C:\Users\ke2n> python -c ‘import obspy; print(obspy.version)’
1.2.2
But if I try
conda create -n ppvmotion python=3 jupyter jupyterlab matplotlib obspy numpy
I get
PackagesNotFoundError: The following packages are not available from current channels:
what I am doing wrong?
Hello kpjamro,
If you have freshly re-installed Anaconda you shouldn’t need this check, but please verify that you are using the latest version (23.3.1). From the command line, you can use conda -V or conda --version to get your current conda version.
After this, you don’t need to import any package in the PPV environment. The command
conda create -n ppvmotion python=3 jupyter jupyterlab matplotlib obspy numpy
will automatically create a new environment named ppvmotion and download/install all required packages/modules. I’ve just tested it on my Anaconda installation on Win11 and everything appears to work fine.
I guess I was trying too hard 
Starting over - it works. Although, I did not figure out where the .ipynb files were.
Instead I just copy/pasted the text of the github file, and filled in the needed spots, and it works on an R3D.
I changed it so that it only asks for one channel - EHZ - and I can make a plot for my system.
Actually, there are three rows in the plot, with just the top graph filled in.
If I set nrows = 1, it throws an error at me. I don’t know how to make only one graph.
But that is pretty good progress for the first day.
Thanks.
OK - total redo of this post after more fiddling.
Maybe Branden or someone at RS can confirm.
I think the way that Obspy should work is that it looks at the instrument XML file finds if the instrument reads in M or M/S or M/S^2 and then integrates or differentiates as needed so that all three outputs can be produced (ACC, VEL or DISP).
yes?
In that case, PPV should simply use the VEL output. The trick with PPV is that, strictly speaking, you need three dimensions and do the Pythagorean sum to find which point in time has the greatest velocity along some diagonal. With an RS1, you only have one axis so it is simple.
The other requirement I have seen is that you say what frequency is associated with the peak. In one reference they suggest finding the waveform zero crossing on either side of the peak and measuring the time between zeros. The inverse of that is the approximate frequency. I don’t know how to do that in Obspy.
BTW - I found why the ACC graph was so low - the sample code had the bandpass filter applied after the conversion to ACC. You need to filter the instrument output before doing the differentiation.
All right kpjamro, lots of progress here, great work!
Let’s proceed step-by-step:
Copy-pasting is more than enough in the end, so it was good thinking on your part.
If you want to plot a single channel (EHZ, in your case) you will have to adapt the code a bit. Just replace the block from line 38 to line 67 with the following:
# generate figure
fig, ax = plt.subplots(figsize=(20, 10), nrows=1, ncols=1)
fig.tight_layout(pad=3.0)
# plot peak particle velocity charts for each channel of the RS3D unit
if len(list(sta)) == 3:
# download and process waveforms for each channel
st = client.get_waveforms(network, station, loc, val, starttime=start_time, endtime=end_time, attach_response=True)
st.merge(method=0, fill_value='latest')
st.detrend(type="demean")
st.filter("bandpass", freqmin=f1, freqmax=f2, corners=cor)
st.remove_response(output="DISP")
# beautify the plots
ax.plot(st[0].times(reftime=start_time), 2*np.pi*freq*st[0].data*1000, lw=1)
ax.margins(x=0.0001)
ax.grid(color='dimgray', ls = '-.', lw = 0.33)
ax.set_ylabel("Peak Particle Velocity (mm/s)")
ax.set_xlabel("Time after Start (s)")
staleg = str(station + " - " + val + " - " + "Filter: bandpass " + str(f1) + " - " + str(f2) + " Hz")
ax.legend([staleg])
# save the final figure
plt.savefig("PPVMotion_" + station + ".png")
else:
# if unit is not a RS3D, warn the user
print("Please insert a RS3D Station.")
The instrument output is always in raw Counts, a unitless number that is essentially an analogue for the amount of voltage on the circuit of the given channel at the given time. (IRIS: Frequently Answered Question ). ObsPy then reads both this data and the instrument response file to output appropriate converted Counts to m (DISP), m/s (VEL), or m/s**2 (ACC).
With the formula that is used in this code from the cited source (https://www.castlegroup.co.uk/ground-vibration/), the a parameter is in mm; thus DISP is required. If, however, you will be using a different formula, feel free to adapt the conversion to suit your needs.
I don’t know how to create this either, but maybe a way could be plotting the spectrogram alongside the waveform and see if, from there, something can be understood. These links could be of use: Plotting Spectrograms — ObsPy 1.4.0 documentation - RPiShakeCode/Seismogram from Truro School start of data logging activity for spectogram.ipynb at master · wmvanstone/RPiShakeCode · GitHub
Thank you for the final comment! This feedback is always useful!
Thanks.
With regard to point 4 the formula: the cited source imagines that your wave is sinusoidal and that you know the frequency of it. With 100 sps the maximum the wave frequency could be is just under 50 Hz. But, in the example, it is filtered down to a maximum of 2 Hz (and may be as low as 0.7). If you plug in 100 Hz for the frequency, you overestimate the velocity by a factor of at least 50 times (100/2) . I won’t comment more on that other than to suggest using “f2” instead of “freq” as a more realistic value for that approach.
In my case, it seemed best just to ask Ospby to give me the velocity directly and use the trace max function to find the maximum value. I produced three graphs - one for displacement, one for velocity (with max - which is the PPV) and one for acceleration (with max):
I would like fewer decimal digits but this will do for now.
I will also scale velocity to be in inches/second which is what local codes require and to express acceleration in fractional percent of gravity - since that seems often used when describing human sensitivity to vibration.
Those are pretty big vibrations for a Raspberry Shake …
Regarding the decimal values, you can round the variable you want to display (PPV, max acceleration, etc.) with this:
round(PPV,2)
where PPV is your variable, and the number after the comma is the number of decimal values you want to display. You can adjust this to suit your needs.
Regarding the rest, this is great work! As it could be of interest to other users too, when you finish optimizing your code, could you post it here in full? So that there will be a good reference point for other interested parties in the future.
Thank you!
Using the round() function gave me an error message
mymax = round(st1.max(),2)
type list doesn’t define round method
np.round is accepted, but it seems that I had to make sure the variable was a single value not a list.
max1 = st1.max()
staleg1 = str(“PPVz " + str(np.round(max1,9)) + " m/sec”)
this actually gives 3 digits when the number is on the order of 1E-7.
For example “6.53e-7”
Hello,
it is possible that the previous mymax value was interpreted as a list for some reason. I cannot be more precise without seeing the code that generates it.
The second part seems to do what you had in mind. I imagine that the PPV values will always be small numbers, so they should match what you need?
putting the np.round() function inside the string building line is what caused the problem.
staleg1 = str(“PPVz " + str(np.round(st1.max(),2)) + " m/sec”)
in principle, st1.max() is just a single number … although I suppose you could have more than one if they were all identical values.
Assigning it to “max1” perhaps picks the first value from the list and then things work.