marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

Merge branch 'feature/magnitude4QtPyLoT'

Browse Source 
Conflicts:
	pylot/core/util/utils.py
This commit is contained in:
Sebastian Wehling-Benatelli
2016-09-13 12:02:21 +02:00
parent 49b07163ad 15700b074d
commit d7680f8ae3
8 changed files with 341 additions and 166 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
pylot/core/analysis/magnitude.py
View File

@@ -14,6 +14,7 @@ from pylot.core.util.utils import getPatternLine
from scipy.optimize import curve_fit
from scipy import integrate, signal
from pylot.core.io.data import Data
from pylot.core.util.dataprocessing import restitute_data
class Magnitude(object):
@@ -304,7 +305,7 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
return Mo, Mw
def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp, iplot):
def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp, iplot):
'''
Subfunction to calculate the source spectrum and to derive from that the plateau
(usually called omega0) and the corner frequency assuming Aki's omega-square
@@ -342,11 +343,8 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
print ("Calculating source spectrum ....")
# get Q value
qu = Qp.split('f**')
# constant Q
Q = int(qu[0])
# A, i.e. power of frequency
A = float(qu[1])
Q, A = qp
delta = delta * 1000 # hypocentral distance in [m]
fc = None
@@ -354,7 +352,7 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
data = Data()
wf_copy = wfstream.copy()
[cordat, restflag] = data.restituteWFData(inventory, wf_copy)
[cordat, restflag] = restitute_data(wf_copy, inventory)
if restflag == 1:
zdat = cordat.select(component="Z")
if len(zdat) == 0:

139
pylot/core/io/data.py
View File

@@ -2,15 +2,13 @@
# -*- coding: utf-8 -*-
import copy
import glob
import os
from obspy import read_events, read_inventory
from obspy import read_events
from obspy.core import read, Stream, UTCDateTime
from obspy.core.event import Event
from obspy.io.xseed import Parser
from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \
picksdict_from_pilot
picksdict_from_pilot, merge_picks
from pylot.core.util.errors import FormatError, OverwriteError
from pylot.core.util.utils import fnConstructor, getGlobalTimes
@@ -269,139 +267,6 @@ class Data(object):
"""
self.get_evt_data().picks = []
def restituteWFData(self, invdlpath, streams=None):
"""
:param invdlpath:
:param streams:
:return:
"""
restflag = 0
if streams is None:
st_raw = self.getWFData()
st = st_raw.copy()
else:
st = streams
for tr in st:
# remove underscores
if len(tr.stats.station) > 3 and tr.stats.station[3] == '_':
tr.stats.station = tr.stats.station[0:3]
dlp = '%s/*.dless' % invdlpath
invp = '%s/*.xml' % invdlpath
respp = '%s/*.resp' % invdlpath
dlfile = glob.glob(dlp)
invfile = glob.glob(invp)
respfile = glob.glob(respp)
# check for dataless-SEED file
if len(dlfile) >= 1:
print("Found dataless-SEED file(s)!")
print("Reading meta data information ...")
for j in range(len(dlfile)):
print("Found dataless-SEED file %s" % dlfile[j])
parser = Parser('%s' % dlfile[j])
for i in range(len(st)):
# check, whether this trace has already been corrected
try:
st[i].stats.processing
except:
try:
print(
"Correcting %s, %s for instrument response "
"..." % (st[i].stats.station,
st[i].stats.channel))
# get corner frequencies for pre-filtering
fny = st[i].stats.sampling_rate / 2
fc21 = fny - (fny * 0.05)
fc22 = fny - (fny * 0.02)
prefilt = [0.5, 0.9, fc21, fc22]
# instrument correction
st[i].simulate(pre_filt=prefilt,
seedresp={'filename': parser,
'date': st[
i].stats.starttime,
'units': "VEL"})
restflag = 1
except ValueError as e:
vmsg = '{0}'.format(e)
print(vmsg)
else:
print("Trace has already been corrected!")
# check for inventory-xml file
if len(invfile) >= 1:
print("Found inventory-xml file(s)!")
print("Reading meta data information ...")
for j in range(len(invfile)):
print("Found inventory-xml file %s" % invfile[j])
inv = read_inventory(invfile[j], format="STATIONXML")
for i in range(len(st)):
# check, whether this trace has already been corrected
try:
st[i].stats.processing
except:
try:
print("Correcting %s, %s for instrument response "
"..." % (st[i].stats.station,
st[i].stats.channel))
# get corner frequencies for pre-filtering
fny = st[i].stats.sampling_rate / 2
fc21 = fny - (fny * 0.05)
fc22 = fny - (fny * 0.02)
prefilt = [0.5, 0.9, fc21, fc22]
# instrument correction
st[i].attach_response(inv)
st[i].remove_response(output='VEL',
pre_filt=prefilt)
restflag = 1
except ValueError as e:
vmsg = '{0}'.format(e)
print(vmsg)
else:
print("Trace has already been corrected!")
# check for RESP-file
if len(respfile) >= 1:
print("Found response file(s)!")
print("Reading meta data information ...")
for j in range(len(respfile)):
print("Found RESP-file %s" % respfile[j])
for i in range(len(st)):
# check, whether this trace has already been corrected
try:
st[i].stats.processing
except:
try:
print("Correcting %s, %s for instrument response "
"..." % (st[i].stats.station,
st[i].stats.channel))
# get corner frequencies for pre-filtering
fny = st[i].stats.sampling_rate / 2
fc21 = fny - (fny * 0.05)
fc22 = fny - (fny * 0.02)
prefilt = [0.5, 0.9, fc21, fc22]
# instrument correction
seedresp = {'filename': respfile[0],
'date': st[0].stats.starttime,
'units': "VEL"}
st[i].simulate(paz_remove=None, pre_filt=prefilt,
seedresp=seedresp)
restflag = 1
except ValueError as e:
vmsg = '{0}'.format(e)
print(vmsg)
else:
print("Trace has already been corrected!")
if len(respfile) < 1 and len(invfile) < 1 and len(dlfile) < 1:
print("No dataless-SEED file,inventory-xml file nor RESP-file "
"found!")
print("Go on processing data without source parameter "
"determination!")
return st, restflag
def get_evt_data(self):
"""

23
pylot/core/io/phases.py
View File

@@ -527,3 +527,26 @@ def writephases(arrivals, fformat, filename):
Ao))
fid.close()
def merge_picks(event, picks):
"""
takes an event object and a list of picks and searches for matching
entries by comparing station name and phase_hint and overwrites the time
and time_errors value of the event picks' with those from the picks
without changing the resource identifiers
:param event: `obspy.core.event.Event` object (e.g. from NLLoc output)
:param picks: list of `obspy.core.event.Pick` objects containing the
original time and time_errors values
:return: merged `obspy.core.event.Event` object
"""
for pick in picks:
time = pick.time
err = pick.time_errors
phase = pick.phase_hint
station = pick.waveform_id.station_code
for p in event.picks:
if p.waveform_id.station_code == station and p.phase_hint == phase:
p.time, p.time_errors = time, err
del time, err, phase, station
return event

10
pylot/core/pick/autopick.py
View File

@@ -17,6 +17,7 @@ from pylot.core.pick.charfuns import CharacteristicFunction
from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
getSNR, fmpicker, checkPonsets, wadaticheck
from pylot.core.util.dataprocessing import restitute_data
from pylot.core.util.utils import getPatternLine
from pylot.core.io.data import Data
from pylot.core.analysis.magnitude import WApp
@@ -564,10 +565,10 @@ def autopickstation(wfstream, pickparam, verbose=False):
hdat = edat.copy()
hdat += ndat
h_copy = hdat.copy()
[cordat, restflag] = data.restituteWFData(invdir, h_copy)
[cordat, restflag] = restitute_data(h_copy, invdir)
# calculate WA-peak-to-peak amplitude
# using subclass WApp of superclass Magnitude
if restflag == 1:
if restflag:
if Sweight < 4:
wapp = WApp(cordat, mpickS, mpickP + sstop, iplot)
else:
@@ -577,6 +578,9 @@ def autopickstation(wfstream, pickparam, verbose=False):
(0.5 * (mpickP + sstop)), iplot)
Ao = wapp.getwapp()
else:
print("Go on processing data without source parameter "
"determination!")
else:
msg = 'Bad initial (AIC) S-pick, skipping this onset!\n' \
@@ -598,7 +602,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
hdat = edat.copy()
hdat += ndat
h_copy = hdat.copy()
[cordat, restflag] = data.restituteWFData(invdir, h_copy)
[cordat, restflag] = restitute_data(h_copy, invdir)
if restflag == 1:
# calculate WA-peak-to-peak amplitude
# using subclass WApp of superclass Magnitude

132
pylot/core/util/dataprocessing.py
View File

@@ -3,9 +3,16 @@
import os
import glob
from obspy import UTCDateTime
import sys
import numpy as np
from obspy import UTCDateTime, read_inventory, read
from obspy.io.xseed import Parser
from pylot.core.util.utils import key_for_set_value, find_in_list, \
remove_underscores
def time_from_header(header):
"""
Function takes in the second line from a .gse file and takes out the date and time from that line.
@@ -148,6 +155,129 @@ def evt_head_check(root_dir, out_dir = None):
print(nfiles)
def restitute_data(data, path_to_inventory, unit='VEL', force=False):
"""
takes a data stream and a path_to_inventory and returns the corrected
waveform data stream
:param data: seismic data stream
:param path_to_inventory: path to inventory folder or file
:param unit: unit to correct for (default: 'VEL')
:param force: force restitution for already corrected traces (default:
False)
:return: corrected data stream
"""
restflag = list()
data = remove_underscores(data)
dlfile = list()
invfile = list()
respfile = list()
inv = dict(dless=dlfile, xml=invfile, resp=respfile)
if os.path.isfile(path_to_inventory):
ext = os.path.splitext(path_to_inventory)[1].split('.')[1]
inv[ext] += [path_to_inventory]
else:
for ext in inv.keys():
inv[ext] += glob.glob1(path_to_inventory, '*.{0}'.format(ext))
invtype = key_for_set_value(inv)
if invtype is None:
print("Neither dataless-SEED file,inventory-xml file nor RESP-file "
"found!")
return data, False
# loop over traces
for tr in data:
seed_id = tr.get_id()
# check, whether this trace has already been corrected
if 'processing' in tr.stats.keys() and not force:
print("Trace {0} has already been corrected!".format(seed_id))
continue
stime = tr.stats.starttime
seedresp = None
prefilt = get_prefilt(tr)
if invtype == 'resp':
fresp = find_in_list(inv[invtype], seed_id)
if not fresp:
raise IOError('no response file found '
'for trace {0}'.format(seed_id))
fname = fresp
seedresp = dict(filename=fname,
date=stime,
units=unit)
kwargs = dict(paz_remove=None, pre_filt=prefilt, seedresp=seedresp)
elif invtype == 'dless':
if len(inv[invtype]) > 1:
fname = Parser(find_in_list(inv[invtype], seed_id))
else:
fname = Parser(inv[invtype][0])
seedresp = dict(filename=fname,
date=stime,
units=unit)
kwargs = dict(pre_filt=prefilt, seedresp=seedresp)
elif invtype == 'xml':
invlist = inv[invtype]
if len(invlist) > 1:
finv = find_in_list(invlist, seed_id)
else:
finv = invlist[0]
inventory = read_inventory(finv, format='STATIONXML')
# apply restitution to data
if invtype in ['resp', 'dless']:
tr.simulate(**kwargs)
else:
tr.attach_response(inventory)
tr.remove_response(output=unit,
pre_filt=prefilt)
restflag.append(True)
# check if ALL traces could be restituted, take care of large datasets
# better try restitution for smaller subsets of data (e.g. station by
# station)
if len(restflag) > 0:
restflag = np.all(restflag)
else:
restflag = False
return data, restflag
def get_prefilt(trace, tlow=(0.5, 0.9), thi=(5., 2.), verbosity=0):
"""
takes a `obspy.core.stream.Trace` object, taper parameters tlow and thi and
returns the pre-filtering corner frequencies for the cosine taper for
further processing
:param trace: seismic data trace
:type trace: `obspy.core.stream.Trace`
:param tlow: tuple or list containing the desired lower corner
frequenices for a cosine taper
:type tlow: tuple or list
:param thi: tuple or list containing the percentage values of the
Nyquist frequency for the desired upper corner frequencies of the cosine
taper
:type thi: tuple or list
:param verbosity: verbosity level
:type verbosity: int
:return: pre-filt cosine taper corner frequencies
:rtype: tuple
..example::
>>> st = read()
>>> get_prefilt(st[0])
(0.5, 0.9, 47.5, 49.0)
"""
if verbosity:
print("Calculating pre-filter values for %s, %s ..." % (
trace.stats.station, trace.stats.channel))
# get corner frequencies for pre-filtering
fny = trace.stats.sampling_rate / 2
fc21 = fny - (fny * thi[0]/100.)
fc22 = fny - (fny * thi[1]/100.)
return (tlow[0], tlow[1], fc21, fc22)
if __name__ == "__main__":
import doctest

31
pylot/core/util/utils.py
View File

@@ -328,6 +328,23 @@ def isIterable(obj):
return False
return True
def key_for_set_value(d):
"""
takes a dictionary and returns the first key for which's value the
boolean is True
:param d: dictionary containing values
:type d: dict
:return: key to the first non-False value found; None if no value's
boolean equals True
"""
r = None
for k, v in d.items():
if v:
return k
return r
def prepTimeAxis(stime, trace):
'''
takes a starttime and a trace object and returns a valid time axis for
@@ -377,6 +394,20 @@ def find_horizontals(data):
return rval
def remove_underscores(data):
"""
takes a `obspy.core.stream.Stream` object and removes all underscores
from stationnames
:param data: stream of seismic data
:type data: `obspy.core.stream.Stream`
:return: data stream
"""
for tr in data:
# remove underscores
tr.stats.station = tr.stats.station.strip('_')
return data
def scaleWFData(data, factor=None, components='all'):
"""
produce scaled waveforms from given waveform data and a scaling factor,