PPSD a spektrální charakteristiky

PPSD = probabilistic power spectral density

In [1]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8

Populating the interactive namespace from numpy and matplotlib

Jak spočítat PPSD

Načti funkce a moduly:

In [2]:
from obspy.signal import PPSD
from obspy import UTCDateTime, read_inventory, read
from obspy.imaging.cm import pqlx
import os

Parametry:

In [3]:
fold_seed = 'SEEDS'
fold_ppsd = 'PPSD_ALPARRAY_02'
period_lim = (0.04, 200)

Staniční metadata:

In [4]:
inv = read_inventory(os.path.join(fold_seed, 'Z3.A071A.xml'))
print(inv.networks[0].stations[0].code, inv.networks[0].stations[0].site.name)
tr = read(os.path.join(fold_seed, 'A071A151030000000.HHE'))[0]

A071A Stare Sedliste, Czech Republic

Inicializace PPSD:

In [5]:
ppsd = PPSD(tr.stats, metadata=inv)
print(len(ppsd._times_processed))

0

Načti signály vstupující do PPSD:

In [6]:
st = read(os.path.join(fold_seed, 'A071A1510*'))
print(st)

2 Trace(s) in Stream:
Z3.A071A..HHE | 2015-10-30T00:00:00.000000Z - 2015-10-30T23:59:59.990000Z | 100.0 Hz, 8640000 samples
Z3.A071A..HHE | 2015-10-31T00:00:00.000000Z - 2015-10-31T23:59:59.990000Z | 100.0 Hz, 8640000 samples

Výpočet PPSD, včetně odstranění přístrojové odezvy:

In [7]:
ppsd.add(st)
print(len(ppsd._times_processed))

94
In [8]:
plt.rcParams['figure.figsize'] = 12, 8
ppsd.plot()

Uložení PPSD:

In [9]:
name = 'ppsd_Z3.A071A..HHE.npz'
ppsd.save_npz(name)
In [10]:
ppsd = []  # vynulujeme

Přišla nová data, joj, už se nám pomíchala se starýma, vše načteme a zapracujeme do PPSD:

In [11]:
st = read(os.path.join(fold_seed, 'A071A15*'))
print(st)

4 Trace(s) in Stream:
Z3.A071A..HHE | 2015-10-30T00:00:00.000000Z - 2015-10-30T23:59:59.990000Z | 100.0 Hz, 8640000 samples
Z3.A071A..HHE | 2015-10-31T00:00:00.000000Z - 2015-10-31T23:59:59.990000Z | 100.0 Hz, 8640000 samples
Z3.A071A..HHE | 2015-11-01T00:00:00.000000Z - 2015-11-01T23:59:59.990000Z | 100.0 Hz, 8640000 samples
Z3.A071A..HHE | 2015-11-02T00:00:00.000000Z - 2015-11-02T23:59:59.990000Z | 100.0 Hz, 8640000 samples

Načtení původně zpracovaných dat:

In [12]:
ppsd = PPSD.load_npz(name, metadata=inv)
print(len(ppsd._times_processed))

94

Přidáme nové, program si sám pohlídá, jestli tam už dané signály jsou (vypíše varovnou hlášku a nevezme je) anebo nejsou (zapracuje je).

Jak se zbavit varovných hlášek:

import warnings warnings.filterwarnings('ignore', 'Already covered time')

In [13]:
ppsd.add(st)

C:\ProgramData\Miniconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py:726: UserWarning: Already covered time spans detected (e.g. 2015-10-30T00:00:00.000000Z), skipping these slices.
  warnings.warn(msg)
C:\ProgramData\Miniconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py:726: UserWarning: Already covered time spans detected (e.g. 2015-10-30T00:30:00.000000Z), skipping these slices.
  warnings.warn(msg)
C:\ProgramData\Miniconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py:726: UserWarning: Already covered time spans detected (e.g. 2015-10-30T01:00:00.000000Z), skipping these slices.
  warnings.warn(msg)
C:\ProgramData\Miniconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py:726: UserWarning: Already covered time spans detected (e.g. 2015-10-30T01:30:00.000000Z), skipping these slices.
  warnings.warn(msg)
...
Out[13]:
True
In [14]:
print(len(ppsd._times_processed))
ppsd.plot()

190

Další možnosti vykreslování

calculate_histogram - zúží histogram jen na vymezené období; zpět pomocí stejného příkazu bez parametrů
plot_temporal - časový vývoj jen pro určené periody; umí vymezit období (viz _stack_selection)
plot_spectrogram - časový vývoj všech period, neumožňuje vymezit období (škoda)

Čtyřdenní spektrogram, vidím 3x denně zvonění, 2 lokální eventy?, mikroseismy, dlouhoperiodický šum pravděpodobně v důsledku změn teploty a tlaku.

In [15]:
ppsd.plot_spectrogram(cmap=pqlx)

Diurnální variace

In [16]:
fin = 'ppsd_Z3.A074A..HHE_2016.npz'
ppsd = PPSD.load_npz(os.path.join(fold_ppsd, fin))
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

Amplitudy spekter s periodama 0.06 s, 0.3 s a 30 s za 4 letní dny:

In [17]:
ppsd.plot_temporal([0.06, 0.3, 30], starttime=UTCDateTime('2016-07-31'), endtime=UTCDateTime('2016-08-05'))

Noční (17-03) a denní režim (04-16):

In [18]:
print('NIGHTLY PPSD')
ppsd.calculate_histogram(time_of_weekday=[(-1,17,24),(-1,0,3)])  # přes noc
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

print('DAYLY PPSD')
ppsd.calculate_histogram(time_of_weekday=[(-1,4,16)])  # přes den
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

ppsd.calculate_histogram() # zruší filtrování

NIGHTLY PPSD

DAYLY PPSD

Sezónní variace

In [19]:
fin = 'ppsd_Z3.A086A..HHZ_2016.npz'
ppsd = PPSD.load_npz(os.path.join(fold_ppsd, fin))
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

Amplitudy spekter s periodama 0.1 s, 1 s a 100 s za celý rok (a v nočních hodinách):

In [20]:
ppsd.plot_temporal([0.1, 1, 100], time_of_weekday=[(-1,0,2)], marker='o', linestyle='')

Letní a zimní režim na dlouhých periodách

In [21]:
print('ZIMA 2015/16, 2. polovina')
ppsd.calculate_histogram(endtime=UTCDateTime('2016-05-01'))
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

print('LETO')
ppsd.calculate_histogram(starttime=UTCDateTime('2016-05-01'), endtime=UTCDateTime('2016-09-15'))
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)

print('ZIMA 2016/17, 1. polovina')
ppsd.calculate_histogram(starttime=UTCDateTime('2016-09-15'))
ppsd.plot(cmap=pqlx, period_lim=period_lim, grid=False)
ppsd.calculate_histogram()

ZIMA 2015/16, 2. polovina

LETO

ZIMA 2016/17, 1. polovina