Oprava raw mseedu


In [1]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
import warnings
warnings.filterwarnings("ignore")
Populating the interactive namespace from numpy and matplotlib
In [2]:
from obspy.core import read, UTCDateTime

Raw mseed je výstup z datalogeru GAIA. Střídají se v něm krátké signály z jednotlivých složek, Z N E Z N ... a uložen je po hodinách. Načteme část takového souboru:

In [3]:
st  = read('07282000.00', endtime=UTCDateTime('2015-07-28T20:50:00'))
In [4]:
print(st)
9 Trace(s) in Stream:
CZ.AAE10..HHE | 2015-07-28T20:00:00.130000Z - 2015-07-28T20:22:11.780000Z | 100.0 Hz, 133166 samples
CZ.AAE10..HHE | 2015-07-28T20:22:14.450900Z - 2015-07-28T20:22:18.560900Z | 100.0 Hz, 412 samples
CZ.AAE10..HHE | 2015-07-28T20:22:15.910000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 166410 samples
CZ.AAE10..HHN | 2015-07-28T20:00:02.300000Z - 2015-07-28T20:22:14.030000Z | 100.0 Hz, 133174 samples
CZ.AAE10..HHN | 2015-07-28T20:22:16.700900Z - 2015-07-28T20:22:20.790900Z | 100.0 Hz, 410 samples
CZ.AAE10..HHN | 2015-07-28T20:22:18.140000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 166187 samples
CZ.AAE10..HHZ | 2015-07-28T20:00:03.670000Z - 2015-07-28T20:22:14.780000Z | 100.0 Hz, 133112 samples
CZ.AAE10..HHZ | 2015-07-28T20:22:17.450900Z - 2015-07-28T20:22:21.520900Z | 100.0 Hz, 408 samples
CZ.AAE10..HHZ | 2015-07-28T20:22:18.870000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 166114 samples

Pro větší přehlednost se zaměříme jen na složku Z a necháme si vypsat časové díry a přesahy.

In [5]:
stZ = st.select(channel="*Z")
print(stZ)
stZ.print_gaps()
3 Trace(s) in Stream:
CZ.AAE10..HHZ | 2015-07-28T20:00:03.670000Z - 2015-07-28T20:22:14.780000Z | 100.0 Hz, 133112 samples
CZ.AAE10..HHZ | 2015-07-28T20:22:17.450900Z - 2015-07-28T20:22:21.520900Z | 100.0 Hz, 408 samples
CZ.AAE10..HHZ | 2015-07-28T20:22:18.870000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 166114 samples
Source            Last Sample                 Next Sample                 Delta           Samples 
CZ.AAE10..HHZ     2015-07-28T20:22:14.780000Z 2015-07-28T20:22:17.450900Z 2.660900        266     
CZ.AAE10..HHZ     2015-07-28T20:22:21.520900Z 2015-07-28T20:22:18.870000Z -2.660900       -266    
Total: 1 gap(s) and 1 overlap(s)
In [6]:
t1 = stZ[1].stats.starttime - 10
t2 = stZ[1].stats.endtime + 10
stZ.plot(starttime=t1, endtime=t2)

Druhá Trace je odskočena, upravíme její počátek tak, aby se zařadila zpět. Pokud se vše povede, příkaz merge(method=-1) udělá ze tří Trace oddělených mezerou či přesahem jednu souvislou stopu.

In [7]:
stZ[1].stats.starttime -= 2.6609
stZ.merge(method=-1)
print(stZ)
stZ.plot(starttime=t1, endtime=t2)
1 Trace(s) in Stream:
CZ.AAE10..HHZ | 2015-07-28T20:00:03.670000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 299634 samples

Tím je zároveň opravena stopa Z i v původním st. Podobně postupujeme i u zbylých složek N, E.

In [8]:
st.select(channel="*N")[1].stats.starttime -= 2.6609
st.select(channel="*E")[1].stats.starttime -= 2.6609
st.merge(method=-1)
print(st)
3 Trace(s) in Stream:
CZ.AAE10..HHZ | 2015-07-28T20:00:03.670000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 299634 samples
CZ.AAE10..HHN | 2015-07-28T20:00:02.300000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 299771 samples
CZ.AAE10..HHE | 2015-07-28T20:00:00.130000Z - 2015-07-28T20:50:00.000000Z | 100.0 Hz, 299988 samples

Ajta krajta, ještě je špatně network.

In [9]:
for tr in st:
    tr.stats.network = 'XT'

Nyní můžeme opravený mseed uložit, pro vytoření názvu souboru použijeme metadata v objektu Trace.

In [10]:
for tr in st:
    output = tr.id + tr.stats.starttime.strftime('.%Y.%j.%H%M%S') + '.seed'
    print(output)
    tr.write(output, format='MSEED')
XT.AAE10..HHZ.2015.209.200003.seed
XT.AAE10..HHN.2015.209.200002.seed
XT.AAE10..HHE.2015.209.200000.seed