Noise

Lab: Probabilistic Power Spectral Densities

Seismo-Live: http://seismo-live.org

- Tobias Megies (@megies)

In [1]:

```
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use("bmh")
plt.rcParams['figure.figsize'] = 10, 6
```

- read waveform data from file
`data/GR.FUR..BHN.D.2015.361`

(station`FUR`

, LMU geophysical observatory in FÃ¼rstenfeldbruck) - read corresponding station metadata from file
`data/station_FUR.stationxml`

- print info on both waveforms and station metadata

In [2]:

```
from obspy import read, read_inventory
st = read("data/GR.FUR..BHN.D.2015.361")
inv = read_inventory("data/station_FUR.stationxml")
print(st)
print(inv)
inv.plot(projection="ortho");
```

- compute probabilistic power spectral densities using
`PPSD`

class from obspy.signal, see http://docs.obspy.org/tutorial/code_snippets/probabilistic_power_spectral_density.html (but use the inventory you read from StationXML as metadata) - plot the processed
`PPSD`

(`plot()`

method attached to`PPSD`

object)

In [3]:

```
from obspy.signal import PPSD
tr = st[0]
ppsd = PPSD(stats=tr.stats, metadata=inv)
ppsd.add(tr)
ppsd.plot()
```

Since longer term stacks would need too much waveform data and take way too long to compute, we prepared one year continuous data preprocessed for a single channel of station `FUR`

to play with..

- load long term pre-computed PPSD from file
`PPSD_FUR_HHN.npz`

using`PPSD`

's`load_npz()`

staticmethod (i.e. it is called directly from the class, not an instance object of the class) - plot the PPSD (default is full time-range, depending on how much data and spread is in the data, adjust
`max_percentage`

option of`plot()`

option) (might take a couple of minutes..!) - do a cumulative plot (which is good to judge non-exceedance percentage dB thresholds)

In [4]:

```
from obspy.signal import PPSD
ppsd = PPSD.load_npz("data/PPSD_FUR_HHN.npz")
```

In [5]:

```
ppsd.plot(max_percentage=10)
ppsd.plot(cumulative=True)
```

- do different stacks of the data using the
`calculate_histogram()`

(see docs!) method of`PPSD`

and visualize them - compare differences in different frequency bands qualitatively (anthropogenic vs. "natural" noise)..
- nighttime stack, daytime stack

- advanced exercise: Use the
`callback`

option and use some crazy custom callback function in`calculate_histogram()`

, e.g. stack together all data from birthdays in your family.. or all German holidays + Sundays in the time span.. or from dates of some bands' concerts on a tour.. etc.

In [6]:

```
ppsd.calculate_histogram(time_of_weekday=[(-1, 0, 2), (-1, 22, 24)])
ppsd.plot(max_percentage=10)
ppsd.calculate_histogram(time_of_weekday=[(-1, 8, 16)])
ppsd.plot(max_percentage=10)
```

- do different stacks of the data using the
`calculate_histogram()`

(see docs!) method of`PPSD`

and visualize them - compare differences in different frequency bands qualitatively (anthropogenic vs. "natural" noise)..
- weekdays stack, weekend stack

In [7]:

```
ppsd.calculate_histogram(time_of_weekday=[(1, 0, 24), (2, 0, 24), (3, 0, 24), (4, 0, 24), (5, 0, 24)])
ppsd.plot(max_percentage=10)
ppsd.calculate_histogram(time_of_weekday=[(6, 0, 24), (7, 0, 24)])
ppsd.plot(max_percentage=10)
```

- do different stacks of the data using the
`calculate_histogram()`

(see docs!) method of`PPSD`

and visualize them - compare differences in different frequency bands qualitatively (anthropogenic vs. "natural" noise)..
- seasonal stacks (e.g. northern hemisphere autumn vs. spring/summer, ...)

In [8]:

```
ppsd.calculate_histogram(month=[10, 11, 12, 1])
ppsd.plot(max_percentage=10)
ppsd.calculate_histogram(month=[4, 5, 6, 7])
ppsd.plot(max_percentage=10)
```

- do different stacks of the data using the
`calculate_histogram()`

(see docs!) method of`PPSD`

and visualize them - compare differences in different frequency bands qualitatively (anthropogenic vs. "natural" noise)..
- stacks by specific month
- maybe even combine several of above restrictions.. (e.g. only nighttime on weekends)

In [ ]:

```
```