Integrating velocity over time for a stream

Usage Questions
1

I’m very new to obspy as my expertise is in marine chemistry. I am working with a high school student to test their idea that traffic volume (cars/time) would be related to seismic noise as measured by a local SchoolShake network.

I have, with help from colleagues, been able to download the stream and extract the velocity for the 15 minute time period where cars passing were counted. I have applied a bandpass to extract the velocity in the 8-12 Hz frequency that would likely result from street traffic. What I think I would like to do is to be able to integrate or sum up all of the velocity to compare to traffic/time to test their hypothesis. But I am not sure how to do this. Could someone suggest how this could be done with ObsPy?

2

from obspy import UTCDateTime
from obspy.clients.fdsn import Client
t1=UTCDateTime(“2022-12-01T01:00:00.713Z”)
t2=UTCDateTime(“2022-12-01T01:15:03.143Z”)
client=Client(‘RASPISHAKE’)
st=client.get_waveforms(network=“AM”,station=“RF395”,location=“00”,channel=“EHZ”,starttime=t1,endtime=t2,attach_response=True)
st.plot()

stvel=st.copy() #make a copy of the stream so you don’t remove the response from the original stream
stvel.remove_response(output=“VEL”) # if you add “plot=True” to this function, it will produce a multi-stage figure showing the steps.
stvel.plot() # now will plot the response corrected seismogram with units of velocity. I think m/s, though it may have changed them to nm/s

stvel.filter(“bandpass”, freqmin=8, freqmax=12)
stvel.plot()

3

You could either look at noise levels in frequency domain (see here for some examples) or you could indeed look at the amplitudes in time domain after applying a bandpass (8-12 Hz seems a bit narrow, I would probably do something like 1-20 Hz, but thats personal choice). In time domain a common measure is the I95 level, i.e. you do all preprocessing (response removal, filtering), you take out a time window (e.g. 10 min or 30 min), then demean the short window and then look at what value is not exceeded by 95% of the samples (i.e. take the absolute, then numpy percentile at 95).

4

One important thing is for your preprocessing you want some data buffer at front and end, so that you dont have artifacts of processing in the data you evaluate, i.e. do the response removal from t1 - 300 to t2 + 300 and then only cut out your window of interest at the end.

5

Thanks for the information Tobias. I appreciate the advice about the time window and artifacts as well. Your Jupyter notebook using PPSD is quite comprehensive. It might be more than I need for this student. I really need some simple and easy to explain measure of the noise for the different stations we counted cars passing.

When you say “(i.e. take the absolute, then numpy percentile at 95)”…I can apply numpy commands to the stream to quantify how noisy the trace is? Could you give a simple example in code for me?

6

Tobias, could you walk me through these steps? I am very new at this.

7

Here’s something to give you an idea:

from obspy import read, read_inventory
import matplotlib.pyplot as plt
import numpy as np

st = read("http://examples.obspy.org/BW.BGLD..EH.D.2010.037")
inv = read_inventory("http://examples.obspy.org/dataless.seed.BW_BGLD")

fig, axes = plt.subplots(nrows=5)

tr = st[0]
tr.trim(endtime=tr.stats.starttime + 6*60)
axes[0].plot(tr.data)

tr.remove_response(inv, pre_filt=(0.2, 0.5, 80, 100))
tr.filter("bandpass", freqmin=1, freqmax=20)
axes[1].plot(tr.data)

# the above is tapered at start/end, cut off that part to not introduce bias
tr.trim(tr.stats.starttime + 30, tr.stats.endtime - 30)
axes[2].plot(tr.data)

t1 = tr.stats.starttime
t2 = t1 + 40

noise_levels = []
while t2 < tr.stats.endtime:
    tr_tmp = tr.slice(t1, t2)
    tr_tmp.detrend("demean")
    tr_tmp.data = np.abs(tr_tmp.data)
    i95 = np.percentile(tr_tmp.data, 95)
    noise_levels.append(i95)
    t1 += 20
    t2 = t1 + 40
axes[3].plot(tr_tmp.data)
axes[3].axhline(i95, color="r")

axes[4].plot(noise_levels, color="r")
for ax in axes:
    ax.grid()
plt.show()

Figure_1
Figure_11459×885 135 KB

8

Thank you Tobias. I have been able to use this code as well as approaching it in a slightly different way (I didn’t know to demean) to arrive at roughly the same answer. What is interesting is that it looks like the 95% in the frequency range for traffic is a little lower the greater our car counts which is a little counterintuitive. It suggests that other factors (people in building, air handling, local geology?) are controlling the noise in the schools at the time of our measurements. A negative result but interesting nonetheless. I think next time we would do the traffic estimates on the weekend when building occupancy and services (HVAC etc) were at a minimum.

Rough data from my analysis without demean:

Station Cars/time Velocity (95%)
SMUS 201 3.51E-07
GNS 122.5 7.49E-06
Hakai 246.5 1.50E-07
Nsaanich 39.5 1.20E-06

My code follows:

import numpy as np
import matplotlib as plt
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
t1=UTCDateTime(“2022-12-09T01:08:00.713Z”)
t2=UTCDateTime(“2022-12-09T01:23:03.143Z”)
client=Client(‘RASPISHAKE’)
st=client.get_waveforms(network=“AM”,station=“R01DF”,location=“00”,channel=“EHZ”,starttime=t1,endtime=t2,attach_response=True)
st.plot()

stvel=st.copy()
stvel.remove_response(output=“VEL”)
stvel.plot()

stvel.filter(“bandpass”, freqmin=1, freqmax=20)
stvel.plot()

stvelab = np.absolute(stvel)
p = np.percentile(stvelab, 95)
p