Staniční metadata a přístrojové odezvy¶
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")
In [3]:
from obspy.core import read, UTCDateTime
from obspy import read_inventory
Funkce read_inventory umožňuje načíst staniční metadata uložená v různých formátech do jednotného objektu Inventory.
Podporované formáty: INVENTORYXML (Arclink), RESP, SC3ML (Seiscomp), SEED, XSEED (SEED ascii verze), STATIONTXT, STATIONXML.
Formát XML je ascii, což umožňuje i přímou editaci souborů. Například CZ_XT_Z3.xml byl vytvořen poskládáním všech STATIONXML souborů pro stanice:
- net CZ: české permanetní stanice, metadata uložená v silo:/arch/service/gfu/dataless/XML_2018MAY
- net XT/CZ: stanice MOBNETu v AlpArray-EASI, metadata vlastní výroby
- net Z3/CZ: stanice MOBNETu v AlpArray, metadata vlastní výroby
In [4]:
inv = read_inventory("CZ_XT_Z3.xml")
print(inv)
In [5]:
inv.plot(projection="ortho", label=False, color_per_network=True),
Out[5]:
In [6]:
plt.rcParams['figure.figsize'] = 12, 8
inv.plot(projection='local', color_per_network=True),
Out[6]:
Příklad použití metadat (lat, long, elev) uložených v Catalog pro výpočet arrayové odezvy / Array response function, pomocí procedury array_transff_wavenumber.
In [7]:
import numpy as np
from obspy.imaging.cm import viridis_r as cmap
from obspy.signal.array_analysis import array_transff_wavenumber
# array coordinates
coords = []
for net in inv.networks:
for stat in net.stations:
coords.append([stat.longitude, stat.latitude, stat.elevation/1000.])
coords = np.asarray(coords)
coords_CZ = []
for net in inv.select(network="CZ").networks:
for stat in net.stations:
coords_CZ.append([stat.longitude, stat.latitude, stat.elevation/1000.])
coords_CZ = np.asarray(coords_CZ)
# set limits for wavenumber differences to analyze
klim = 0.1; kstep = klim / 100.
kxmin = -klim; kxmax = klim
kymin = -klim; kymax = klim
# compute transfer function as a function of wavenumber difference
transff = array_transff_wavenumber(coords, klim, kstep, coordsys='lonlat')
transff_CZ = array_transff_wavenumber(coords_CZ, klim, kstep, coordsys='lonlat')
# FIGURE
fig = plt.figure(figsize=(14.5,6))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
# plot CZ + AlpArray
sc1 = ax1.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
transff.T, cmap=cmap)
plt.colorbar(sc1, ax=ax1)
sc1.set_clim(vmin=0., vmax=1.)
ax1.set_xlim(kxmin, kxmax); ax1.set_ylim(kymin, kymax)
ax1.set_xlabel('$K=1/\lambda \quad [m^{-1}]$')
ax1.set_ylabel('$K=1/\lambda \quad [m^{-1}]$')
ax1.set_title('Array response function: CZ + AlpArray')
# plot CZ PERM only
sc2 = plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
transff_CZ.T, cmap=cmap)
plt.colorbar(sc2, ax=ax2)
sc2.set_clim(vmin=0., vmax=1.)
ax2.set_xlim(kxmin, kxmax); plt.ylim(kymin, kymax)
ax2.set_xlabel('$K=1/\lambda \quad [m^{-1}]$')
ax2.set_ylabel('$K=1/\lambda \quad [m^{-1}]$')
ax2.set_title('Array response function - CZ permanent')
plt.tight_layout()
plt.show()
Odezva seismometru¶
Informace k odezvě seismometru a datalogeru jsou uložené v Channel.
Podle metadat má stanice KHC v současnosti tyto instrumenty:
- Streckeisen STS-2 + Quanterra Q330S (BH, HH)
- krátkoperiodický Vegik + Earth Data PR6-24 (SH) / Quanterra Q330S (EH)
In [8]:
inv_KHC_Z = inv.select(station="KHC", channel="*Z")
print(inv_KHC_Z.select(time=UTCDateTime()).networks[0].stations[0].channels[0])
print(inv_KHC_Z.select(time=UTCDateTime()).networks[0].stations[0].channels[0].response)
print(inv_KHC_Z.select(time=UTCDateTime()).networks[0].stations[0].channels[0].response.response_stages[0])
Jeden obrázek lepší než plno textu ...
In [9]:
# díky time=UTCDateTime() vybere jen současné odezvy
inv_KHC_Z.plot_response(min_freq=0.001, time=UTCDateTime(), unwrap_phase=True),
Out[9]:
Pro srovnání stanice GEC2, v současnosti instrumenty:
- Streckeisen STS-2 + Nanometrics HRD24 (BH, HH)
- Geotech GS13 + Nanometrics HRD24 (SH)
In [10]:
inv_GEC2_Z = read_inventory('GR.GEC2.xml', format="STATIONXML").select(channel="*Z")
inv_GEC2_Z.plot_response(min_freq=0.001, time=UTCDateTime()),
Out[10]:
Odstranění odezvy seismometru, simulace jiného seismometru
attach_response- přiřadí ke složce signálu správnou staniční odezvu (už to není povinné)remove_response- odstraní ze signálu staniční odezvusimulateasimulate_seismometer- [odstraní ze signálu staniční odezvu a] přidá jinou odezvu
Během odstraňování odezvy/simulace je důležité, aby filtr byl stabilní, tj. aby se neúměrně nezvýraznily některé frekvence. Toto se dá kontrolovat pomocí další filtrace (pre_filt) a/nebo parametru water_level. Pro lepší přehled ObsPy umožňuje vizualizaci jednotlivých kroků, viz remove_response s parametrem plot=True.
In [11]:
t1 = UTCDateTime('20150806092228'); t2 = t1 + 60*60
st = read('KHC.EHZ.*', starttime=t1, endtime=t2)
pre_filt = [0.001, 0.005, 30, 50]
st.remove_response(inventory=inv, pre_filt=pre_filt, output="VEL", water_level=60, plot=True)
Out[11]:
Příprava signálů na korelační pickování teleseismické P vlny¶
Signály P vlny na všech stanicích musí být podobné, odstranění odezvy seismometru pomůže (obzvláště při kombinaci krátkoperiodických a širokopásmových stanic. Výbornou kontrolou správného postupu jsou stanice s více instrumenty. Zde jako příklad ukážeme KHC a GEC2.
1. Původní signál¶
In [12]:
tP_KHC = UTCDateTime('2015-08-06T09:34:45.97') + 0.86 # manuálně určený začátek P vlny u KHC
tP_GEC2 = UTCDateTime('2015-08-06T09:34:46.84') + 0.52 - 0.076 # manuálně určený začátek P vlny u GEC2
inv_duo = inv.select(station='KHC', channel='*HZ') + inv_GEC2_Z # staniční metadata jen pro tyto stanice
st = read('*.mseed', starttime=tP_KHC-60, endtime=tP_KHC+60).detrend('linear')
print(st)
st.plot(starttime=tP_KHC-5, endtime=tP_KHC+5, equal_scale=False, size=(1200, 600)) # okno +- 5 s kolem P onsetu
2. Simulace WWSSN SP odezvy (bez korekce staničních odezev)¶
aneb bandpass filtr kolem 1 s
In [13]:
paz_wwssn_sp = {'poles': [-3.367788+3.731514j, -3.367788-3.731514j,
-7.037168+4.545562j, -7.037168-4.545562j],
'zeros': [0j, 0j],
'gain': 70.18386,
'sensitivity': 1.}
# sesazení počátků P vln u KHC a GEC2 stanice
for tr in st.select(station='GEC2'):
tr.stats.starttime -= (tP_GEC2 - tP_KHC)
In [14]:
st_simul1 = st.copy()
st_simul1.simulate(paz_remove=None, paz_simulate=paz_wwssn_sp)
st_simul1.plot(starttime=tP_KHC-5, endtime=tP_KHC+5, equal_scale=False, size=(1200, 600))
3. Korekce staničních odezev a simulace WWSSN SP¶
In [15]:
pre_filt = [0.001, 0.005, 2, 4]
st_simul2 = st.copy()
st_simul2.attach_response(inv_duo)
st_simul2.remove_response(pre_filt=pre_filt, water_level=1e6)
st_simul2.simulate(paz_remove=None, paz_simulate=paz_wwssn_sp)
st_simul2.plot(starttime=tP_KHC-5, endtime=tP_KHC+5, equal_scale=False, size=(1200, 600))
Korekce mírných posuntí P vln pomocí korelačního xcorr_pick_correction.
In [16]:
from obspy.signal.cross_correlation import xcorr_pick_correction
# všechny stopy vstupující do xcorr_pick_correction musí mít stejný sampling!
for tr in st_simul1:
if tr.stats.sampling_rate != 100.:
tr.interpolate(sampling_rate=100.)
for tr in st_simul2:
if tr.stats.sampling_rate != 100.:
tr.interpolate(sampling_rate=100.)
In [17]:
#tr1 = st_simul2.select(id='GR.GEC2..HHZ')[0]; tr2 = st_simul2.select(id='GR.GEC2..SHZ')[0]
#tr1 = st_simul2.select(id='GR.GEC2..HHZ')[0]; tr2 = st_simul2.select(id='CZ.KHC..HHZ')[0]
tr1 = st_simul2.select(id='CZ.KHC..HHZ')[0]; tr2 = st_simul2.select(id='CZ.KHC..EHZ')[0];
dt, coeff = xcorr_pick_correction(tP_KHC, tr1, tP_KHC, tr2, 0.5, 2., 1., plot=True)
print('{:s} - {:s} : {:.3f} s {:.0f} %'.format(tr2.id, tr1.id, dt, 100.*coeff))
In [18]:
print('Bez korekce na odezvu seismometru:')
for i,j in [[0,1],[0,3],[3,2]]:
tr1 = st_simul1[i]; tr2 = st_simul1[j]
dt, coeff = xcorr_pick_correction(tP_KHC, tr1, tP_KHC, tr2, 0.5, 2., 1., plot=False)
print('{:<12s} - {:<12s} : {:6.2f} s {:3.0f} %'.format(tr2.id, tr1.id, dt, 100.*coeff))
print('\nPo korekci na odezvu seismometru:')
for i,j in [[0,1],[0,3],[3,2]]:
tr1 = st_simul2[i]; tr2 = st_simul2[j]
dt, coeff = xcorr_pick_correction(tP_KHC, tr1, tP_KHC, tr2, 0.5, 2., 1., plot=False)
print('{:<12s} - {:<12s} : {:6.2f} s {:3.0f} %'.format(tr2.id, tr1.id, dt, 100.*coeff))
ZDE NOTEBOOK UKONČÍME
In [ ]:
# dalsi namety:
# PRU ... SH složka stejně špatná jako u KHC
# A077A ... CMG-3ESP/30s, ukáže rozdíl mezi GAIA1/25Hz a GAIA3
#inv.select(station="A077A", channel="*HZ").plot_response(min_freq=0.001),
In [ ]:
# pomocne funkce k zobrazení response stage
def plot_response_stages(inventory, min_freq, output="VEL", network="*", station="*",
location="*", channel="*", time=None, starttime=None,
endtime=None, axes=None, unwrap_phase=False, show=True,
outfile=None):
import matplotlib.pyplot as plt
if axes:
ax1, ax2 = axes
fig = ax1.figure
else:
fig = plt.figure(figsize=(12,11))
ax1 = plt.subplot2grid((3,1), (0,0), rowspan=2)
ax2 = plt.subplot2grid((3,1), (2,0), sharex=ax1)
matching = inventory.select(network=network, station=station,
location=location, channel=channel, time=time,
starttime=starttime, endtime=endtime)
if len(matching.networks) > 1:
print('WARNING: too much networks, we take the first one from:')
contents = matching.get_contents()
print(" ".join(["%s" % _i for _i in contents["networks"]]))
if len(matching.networks[0].stations) > 1:
print('WARNING: too much stations, we take the first one from:')
contents = matching.networks[0].get_contents()
print(" ".join(["%s" % _i for _i in contents["stations"]]))
if len(matching.networks[0].stations[0].channels) > 1:
print('WARNING: too much channels, we take the first one from:')
print(matching.networks[0].stations[0].__str__().split('\t')[-1])
n_stages = len(matching.networks[0].stations[0].channels[0].response.response_stages)
stat_id = ".".join((matching.networks[0].code, matching.networks[0].stations[0].code,
matching.networks[0].stations[0].channels[0].location_code,
matching.networks[0].stations[0].channels[0].code))
for i in range(1,n_stages+1):
matching.networks[0].stations[0].channels[0].plot(min_freq=min_freq, output=output,
start_stage=i, end_stage=i, axes=(ax1, ax2),
label="Stage {0:d}".format(i), unwrap_phase=unwrap_phase, show=False,
outfile=None)
matching.networks[0].stations[0].channels[0].plot(min_freq=min_freq, output=output,
axes=(ax1, ax2),
label=stat_id, unwrap_phase=unwrap_phase, show=False,
outfile=None)
# final adjustments to plot if we created the figure in here
if not axes:
from obspy.station.response import _adjust_bode_plot_figure
_adjust_bode_plot_figure(fig, show=False)
if outfile:
fig.savefig(outfile)
else:
if show:
plt.show()
print(stat_id)
print(matching.networks[0].stations[0].channels[0].response)
return fig
def inventory_ids(inv):
inv_ids = []
for network in inv.networks:
for station in network.stations:
for channel in station.channels:
if channel.start_date:
start_date = channel.start_date.strftime('%y%m%d%H%M')
else:
start_date = '0000000000'
if channel.end_date:
end_date = channel.end_date.strftime('%y%m%d%H%M')
else:
end_date = '9999999999'
inv_id = ".".join((network.code, station.code, channel.location_code, channel.code,
start_date, end_date))
inv_ids.append(inv_id)
return inv_ids
def inventory_table(inv, output=None):
if output:
f = open(output, 'w')
for network in inv.networks:
for station in network.stations:
for channel in station.channels:
if channel.start_date:
start_date = channel.start_date.strftime('%y%m%d%H%M')
else:
start_date = '0000000000'
if channel.end_date:
end_date = channel.end_date.strftime('%y%m%d%H%M')
else:
end_date = '9999999999'
sss = "{:2s} {:<5s} {:2s} {:8.4f} {:8.4f} {:4.0f} {:3s} {:10s} {:10s} {:3.0f} {:3.0f}".format(
network.code, station.code, channel.location_code,
channel.latitude, channel.longitude, channel.elevation,
channel.code, start_date, end_date, channel.dip, channel.azimuth)
#inv_id += '.{:03d}.{:03d}'.format(int(round(channel.azimuth)), int(round(channel.dip)))
print(sss)
if output:
f.write(sss + '\n')
if output: f.close()
return
def inventory_response_report(inv, output=None, label='---'):
inv_ids = sorted(inventory_ids(inv))
if output:
f = open(output, 'w')
for inv_id in inv_ids:
_inv = inv.select(network=inv_id.split('.')[0], station=inv_id.split('.')[1],
location=inv_id.split('.')[2], channel=inv_id.split('.')[3],
time=UTCDateTime('20'+inv_id.split('.')[4]+'00')+60)
__inv = inv.select(network=inv_id.split('.')[0], station=inv_id.split('.')[1],
location=inv_id.split('.')[2], channel=inv_id.split('.')[3])
sss = '---- '
sss += inv_id
sss += ' -------------------------------------------' + label + '---\n'
sss += _inv.networks[0].stations[0].channels[0].response.__str__()
print(sss)
if output:
f.write(sss + '\n')
if output: f.close()
return