Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

This commit is contained in:
Marcel Paffrath 2016-05-25 14:24:51 +02:00
commit ecf4a7ced7
11 changed files with 474 additions and 294 deletions

View File

@ -51,7 +51,8 @@ from pylot.core.util.errors import FormatError, DatastructureError, \
OverwriteError OverwriteError
from pylot.core.util.connection import checkurl from pylot.core.util.connection import checkurl
from pylot.core.util.utils import fnConstructor, createEvent, getLogin, \ from pylot.core.util.utils import fnConstructor, createEvent, getLogin, \
createCreationInfo, getGlobalTimes getGlobalTimes
from pylot.core.io.location import create_creation_info, create_event
from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \
MPLWidget, PropertiesDlg, HelpForm, createAction, PickDlg MPLWidget, PropertiesDlg, HelpForm, createAction, PickDlg
from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.structure import DATASTRUCTURE
@ -829,8 +830,8 @@ class MainWindow(QMainWindow):
new = NewEventDlg() new = NewEventDlg()
if new.exec_() != QDialog.Rejected: if new.exec_() != QDialog.Rejected:
evtpar = new.getValues() evtpar = new.getValues()
cinfo = createCreationInfo(agency_id=self.agency) cinfo = create_creation_info(agency_id=self.agency)
event = createEvent(evtpar['origintime'], cinfo) event = create_event(evtpar['origintime'], cinfo)
self.data = Data(self, evtdata=event) self.data = Data(self, evtdata=event)
self.setDirty(True) self.setDirty(True)

View File

@ -1 +1 @@
1508-dirty 7c5a-dirty

View File

@ -395,7 +395,7 @@ class SeismicShot(object):
self.picks[traceID]['spe']) = earllatepicker(self.getSingleStream(traceID), self.picks[traceID]['spe']) = earllatepicker(self.getSingleStream(traceID),
nfac, (tnoise, tgap, tsignal), nfac, (tnoise, tgap, tsignal),
self.getPickIncludeRemoved(traceID), self.getPickIncludeRemoved(traceID),
stealthMode=True) stealth_mode=True)
if self.picks[traceID]['epp'] < 0: if self.picks[traceID]['epp'] < 0:
self.picks[traceID]['epp'] self.picks[traceID]['epp']

View File

@ -36,7 +36,7 @@ class AutoPickParameter(object):
========== ========== ======================================= ========== ========== =======================================
''' '''
def __init__(self, fnin=None, fnout=None, **kwargs): def __init__(self, fnin=None, fnout=None, verbosity=0, **kwargs):
''' '''
Initialize parameter object: Initialize parameter object:
@ -60,7 +60,8 @@ class AutoPickParameter(object):
parspl = line.split('\t')[:2] parspl = line.split('\t')[:2]
parFileCont[parspl[0].strip()] = parspl[1] parFileCont[parspl[0].strip()] = parspl[1]
except IndexError as e: except IndexError as e:
self._printParameterError(e) if verbosity > 0:
self._printParameterError(e)
inputFile.seek(0) inputFile.seek(0)
lines = inputFile.readlines() lines = inputFile.readlines()
for line in lines: for line in lines:

222
pylot/core/io/location.py Normal file
View File

@ -0,0 +1,222 @@
from obspy import UTCDateTime
from obspy.core import event as ope
from pylot.core.util.utils import getLogin, getHash
def create_amplitude(pickID, amp, unit, category, cinfo):
'''
:param pickID:
:param amp:
:param unit:
:param category:
:param cinfo:
:return:
'''
amplitude = ope.Amplitude()
amplitude.creation_info = cinfo
amplitude.generic_amplitude = amp
amplitude.unit = ope.AmplitudeUnit(unit)
amplitude.type = ope.AmplitudeCategory(category)
amplitude.pick_id = pickID
return amplitude
def create_arrival(pickresID, cinfo, phase, azimuth=None, dist=None):
'''
create_arrival - function to create an Obspy Arrival
:param pickresID: Resource identifier of the created pick
:type pickresID: :class: `~obspy.core.event.ResourceIdentifier` object
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param phase: name of the arrivals seismic phase
:type phase: str
:param azimuth: azimuth between source and receiver
:type azimuth: float or int, optional
:param dist: distance between source and receiver
:type dist: float or int, optional
:return: An ObsPy :class: `~obspy.core.event.Arrival` object
'''
arrival = ope.Arrival()
arrival.creation_info = cinfo
arrival.pick_id = pickresID
arrival.phase = phase
if azimuth is not None:
arrival.azimuth = float(azimuth) if azimuth > -180 else azimuth + 360.
else:
arrival.azimuth = azimuth
arrival.distance = dist
return arrival
def create_creation_info(agency_id=None, creation_time=None, author=None):
'''
:param agency_id:
:param creation_time:
:param author:
:return:
'''
if author is None:
author = getLogin()
if creation_time is None:
creation_time = UTCDateTime()
return ope.CreationInfo(agency_id=agency_id, author=author,
creation_time=creation_time)
def create_event(origintime, cinfo, originloc=None, etype=None, resID=None,
authority_id=None):
'''
create_event - funtion to create an ObsPy Event
:param origintime: the events origintime
:type origintime: :class: `~obspy.core.utcdatetime.UTCDateTime` object
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param originloc: tuple containing the location of the origin
(LAT, LON, DEP) affiliated with the event which is created
:type originloc: tuple, list
:param etype: Event type str object. converted via ObsPy to a valid event
type string.
:type etype: str
:param resID: Resource identifier of the created event
:type resID: :class: `~obspy.core.event.ResourceIdentifier` object, str
:param authority_id: name of the institution carrying out the processing
:type authority_id: str
:return: An ObsPy :class: `~obspy.core.event.Event` object
'''
etype = ope.EventType(etype)
if originloc is not None:
o = create_origin(origintime, cinfo,
originloc[0], originloc[1], originloc[2])
else:
o = None
if etype is None:
etype = ope.EventType('earthquake') # defaults to 'earthquake'
if not resID:
resID = create_resourceID(origintime, etype, authority_id)
elif isinstance(resID, str):
resID = create_resourceID(origintime, etype, authority_id, resID)
elif not isinstance(resID, ope.ResourceIdentifier):
raise TypeError("unsupported type(resID) for resource identifier "
"generation: %s" % type(resID))
event = ope.Event(resource_id=resID)
event.creation_info = cinfo
event.event_type = etype
if o:
event.origins = [o]
return event
def create_magnitude(originID, cinfo):
'''
create_magnitude - function to create an ObsPy Magnitude object
:param originID:
:type originID:
:param cinfo:
:type cinfo:
:return:
'''
magnitude = ope.Magnitude()
magnitude.creation_info = cinfo
magnitude.origin_id = originID
return magnitude
def create_origin(origintime, cinfo, latitude, longitude, depth):
'''
create_origin - function to create an ObsPy Origin
:param origintime: the origins time of occurence
:type origintime: :class: `~obspy.core.utcdatetime.UTCDateTime` object
:param cinfo:
:type cinfo:
:param latitude: latitude in decimal degree of the origins location
:type latitude: float
:param longitude: longitude in decimal degree of the origins location
:type longitude: float
:param depth: hypocentral depth of the origin
:type depth: float
:return: An ObsPy :class: `~obspy.core.event.Origin` object
'''
assert isinstance(origintime, UTCDateTime), "origintime has to be " \
"a UTCDateTime object, but " \
"actually is of type " \
"'%s'" % type(origintime)
origin = ope.Origin()
origin.time = origintime
origin.creation_info = cinfo
origin.latitude = latitude
origin.longitude = longitude
origin.depth = depth
return origin
def create_pick(origintime, picknum, picktime, eventnum, cinfo, phase, station,
wfseedstr, authority_id):
'''
create_pick - function to create an ObsPy Pick
:param origintime:
:type origintime:
:param picknum: number of the created pick
:type picknum: int
:param picktime:
:type picktime:
:param eventnum: human-readable event identifier
:type eventnum: str
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param phase: name of the arrivals seismic phase
:type phase: str
:param station: name of the station at which the seismic phase has been
picked
:type station: str
:param wfseedstr: A SEED formatted string of the form
network.station.location.channel in order to set a referenced waveform
:type wfseedstr: str, SEED formatted
:param authority_id: name of the institution carrying out the processing
:type authority_id: str
:return: An ObsPy :class: `~obspy.core.event.Pick` object
'''
pickID = eventnum + '_' + station.strip() + '/{0:03d}'.format(picknum)
pickresID = create_resourceID(origintime, 'pick', authority_id, pickID)
pick = ope.Pick()
pick.resource_id = pickresID
pick.time = picktime
pick.creation_info = cinfo
pick.phase_hint = phase
pick.waveform_id = ope.ResourceIdentifier(id=wfseedstr, prefix='file:/')
return pick
def create_resourceID(timetohash, restype, authority_id=None, hrstr=None):
'''
:param timetohash:
:type timetohash
:param restype: type of the resource, e.g. 'orig', 'earthquake' ...
:type restype: str
:param authority_id: name of the institution carrying out the processing
:type authority_id: str, optional
:param hrstr:
:type hrstr:
:return:
'''
assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \
"UTCDateTime object"
hid = getHash(timetohash)
if hrstr is None:
resID = ope.ResourceIdentifier(restype + '/' + hid[0:6])
else:
resID = ope.ResourceIdentifier(restype + '/' + hrstr + '_' + hid[0:6])
if authority_id is not None:
resID.convertIDToQuakeMLURI(authority_id=authority_id)
return resID

View File

@ -1,16 +1,18 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import os
import glob import glob
import warnings
import scipy.io as sio
import obspy.core.event as ope import obspy.core.event as ope
import os
import scipy.io as sio
import warnings
from obspy.core import UTCDateTime from obspy.core import UTCDateTime
from pylot.core.io.location import create_arrival, create_event, \
create_magnitude, create_origin, create_pick
from pylot.core.pick.utils import select_for_phase from pylot.core.pick.utils import select_for_phase
from pylot.core.util.utils import getOwner, createPick, createArrival, \ from pylot.core.util.utils import getOwner, getGlobalTimes
createEvent, createOrigin, createMagnitude, getGlobalTimes
def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
""" """
@ -75,14 +77,14 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
stations = [stat for stat in phases['stat'][0:-1:3]] stations = [stat for stat in phases['stat'][0:-1:3]]
event = createEvent(eventDate, loccinfo, etype='earthquake', resID=eventNum, event = create_event(eventDate, loccinfo, etype='earthquake', resID=eventNum,
authority_id=authority_id) authority_id=authority_id)
lat = float(loc['LAT']) lat = float(loc['LAT'])
lon = float(loc['LON']) lon = float(loc['LON'])
dep = float(loc['DEP']) dep = float(loc['DEP'])
origin = createOrigin(eventDate, loccinfo, lat, lon, dep) origin = create_origin(eventDate, loccinfo, lat, lon, dep)
for n, pick in enumerate(phases['Ptime']): for n, pick in enumerate(phases['Ptime']):
if pick[0] > 0: if pick[0] > 0:
kwargs = {'year': int(pick[0]), kwargs = {'year': int(pick[0]),
@ -115,15 +117,15 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
wffn = os.path.join(sdir, '{0}*{1}*'.format( wffn = os.path.join(sdir, '{0}*{1}*'.format(
stations[n].strip(), '[ne]')) stations[n].strip(), '[ne]'))
print(wffn) print(wffn)
pick = createPick(eventDate, np, picktime, eventNum, pickcinfo, pick = create_pick(eventDate, np, picktime, eventNum, pickcinfo,
phase, stations[n], wffn, authority_id) phase, stations[n], wffn, authority_id)
event.picks.append(pick) event.picks.append(pick)
pickID = pick.get('id') pickID = pick.get('id')
arrival = createArrival(pickID, pickcinfo, phase) arrival = create_arrival(pickID, pickcinfo, phase)
origin.arrivals.append(arrival) origin.arrivals.append(arrival)
np += 1 np += 1
magnitude = createMagnitude(origin.get('id'), loccinfo) magnitude = create_magnitude(origin.get('id'), loccinfo)
magnitude.mag = float(loc['Mnet']) magnitude.mag = float(loc['Mnet'])
magnitude.magnitude_type = 'Ml' magnitude.magnitude_type = 'Ml'
@ -232,7 +234,7 @@ def picksdict_from_picks(evt):
def picks_from_picksdict(picks): def picks_from_picksdict(picks):
picks_list = list() picks_list = list()
for station, onsets in picks.items(): for station, onsets in picks.items():
print('Reading picks on station %s' % station) #print('Reading picks on station %s' % station)
for label, phase in onsets.items(): for label, phase in onsets.items():
if not isinstance(phase, dict) or len(phase) < 3: if not isinstance(phase, dict) or len(phase) < 3:
continue continue
@ -262,21 +264,28 @@ def picks_from_picksdict(picks):
else: else:
pick.polarity = 'undecidable' pick.polarity = 'undecidable'
except KeyError as e: except KeyError as e:
print(e.message, 'No polarity information found for %s' % phase) if 'fm' in e.message: # no polarity information found for this phase
pass
else:
raise e
picks_list.append(pick) picks_list.append(pick)
return picks_list return picks_list
def reassess_pilot_db(root_dir, out_dir=None, fn_param=None): def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0):
import glob import glob
evt_list = glob.glob1(root_dir,'e????.???.??') db_root = os.path.join(root_dir, db_dir)
evt_list = glob.glob1(db_root,'e????.???.??')
for evt in evt_list: for evt in evt_list:
reassess_pilot_event(root_dir, evt, out_dir, fn_param) if verbosity > 0:
print('Reassessing event {0}'.format(evt))
reassess_pilot_event(root_dir, db_dir, evt, out_dir, fn_param, verbosity)
def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbosity=0):
def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None, verbosity=0):
from obspy import read from obspy import read
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import AutoPickParameter
@ -286,15 +295,19 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
import pylot.core.util.defaults as defaults import pylot.core.util.defaults as defaults
fn_param = defaults.AUTOMATIC_DEFAULTS fn_param = defaults.AUTOMATIC_DEFAULTS
default = AutoPickParameter(fn_param) default = AutoPickParameter(fn_param, verbosity)
search_base = os.path.join(root_dir, event_id) search_base = os.path.join(root_dir, db_dir, event_id)
phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat')) phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat'))
if not phases_file: if not phases_file:
return return
print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0])) if verbosity > 1:
print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0]))
picks_dict = picks_from_pilot(phases_file[0]) picks_dict = picks_from_pilot(phases_file[0])
print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict)) if verbosity > 0:
print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict))
datacheck = list()
info = None
for station in picks_dict.keys(): for station in picks_dict.keys():
fn_pattern = os.path.join(search_base, '{0}*'.format(station)) fn_pattern = os.path.join(search_base, '{0}*'.format(station))
try: try:
@ -302,6 +315,21 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
except TypeError as e: except TypeError as e:
print(e.message) print(e.message)
st = read(fn_pattern) st = read(fn_pattern)
except ValueError as e:
if e.message == 'second must be in 0..59':
info = 'A known Error was raised. Please find the list of corrupted files and double-check these files.'
datacheck.append(fn_pattern + ' (time info)\n')
continue
else:
raise ValueError(e.message)
except Exception as e:
if 'No file matching file pattern:' in e.message:
if verbosity > 0:
warnings.warn('no waveform data found for station {station}'.format(station=station), RuntimeWarning)
datacheck.append(fn_pattern + ' (no data)\n')
continue
else:
raise e
for phase in picks_dict[station].keys(): for phase in picks_dict[station].keys():
try: try:
mpp = picks_dict[station][phase]['mpp'] mpp = picks_dict[station][phase]['mpp']
@ -310,10 +338,9 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
continue continue
sel_st = select_for_phase(st, phase) sel_st = select_for_phase(st, phase)
if not sel_st: if not sel_st:
raise warnings.formatwarning( msg = 'no waveform data found for station {station}'.format(station=station)
'no waveform data found for station {station}'.format( warnings.warn(msg, RuntimeWarning)
station=station), category=RuntimeWarning) continue
print(sel_st)
stime, etime = getGlobalTimes(sel_st) stime, etime = getGlobalTimes(sel_st)
rel_pick = mpp - stime rel_pick = mpp - stime
epp, lpp, spe = earllatepicker(sel_st, epp, lpp, spe = earllatepicker(sel_st,
@ -321,7 +348,7 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
default.get('tsnrz' if phase == 'P' else 'tsnrh'), default.get('tsnrz' if phase == 'P' else 'tsnrh'),
Pick1=rel_pick, Pick1=rel_pick,
iplot=None, iplot=None,
) stealth_mode=True)
if epp is None or lpp is None: if epp is None or lpp is None:
continue continue
epp = stime + epp epp = stime + epp
@ -332,13 +359,24 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
if mpp - epp < min_diff: if mpp - epp < min_diff:
epp = mpp - min_diff epp = mpp - min_diff
picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe) picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe)
if datacheck:
if info:
if verbosity > 0:
print(info + ': {0}'.format(search_base))
fncheck = open(os.path.join(search_base, 'datacheck_list'), 'w')
fncheck.writelines(datacheck)
fncheck.close()
del datacheck
# create Event object for export # create Event object for export
evt = ope.Event(resource_id=event_id) evt = ope.Event(resource_id=event_id)
evt.picks = picks_from_picksdict(picks_dict) evt.picks = picks_from_picksdict(picks_dict)
# write phase information to file # write phase information to file
if not out_dir: if not out_dir:
fnout_prefix = os.path.join(root_dir, event_id, '{0}.'.format(event_id)) fnout_prefix = os.path.join(root_dir, db_dir, event_id, '{0}.'.format(event_id))
else: else:
out_dir = os.path.join(out_dir, db_dir)
if not os.path.isdir(out_dir):
os.makedirs(out_dir)
fnout_prefix = os.path.join(out_dir, '{0}.'.format(event_id)) fnout_prefix = os.path.join(out_dir, '{0}.'.format(event_id))
evt.write(fnout_prefix + 'xml', format='QUAKEML') evt.write(fnout_prefix + 'xml', format='QUAKEML')
#evt.write(fnout_prefix + 'cnv', format='VELEST') #evt.write(fnout_prefix + 'cnv', format='VELEST')

View File

@ -16,7 +16,7 @@ import numpy as np
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode=False): def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False):
''' '''
Function to derive earliest and latest possible pick after Diehl & Kissling (2009) Function to derive earliest and latest possible pick after Diehl & Kissling (2009)
as reasonable uncertainties. Latest possible pick is based on noise level, as reasonable uncertainties. Latest possible pick is based on noise level,
@ -45,8 +45,9 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode=False):
LPick = None LPick = None
EPick = None EPick = None
PickError = None PickError = None
if stealthMode is False: if stealth_mode is False:
print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...' print('earllatepicker: Get earliest and latest possible pick'
' relative to most likely pick ...')
x = X[0].data x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
@ -62,8 +63,9 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode=False):
ilup, = np.where(x[isignal] > nlevel) ilup, = np.where(x[isignal] > nlevel)
ildown, = np.where(x[isignal] < -nlevel) ildown, = np.where(x[isignal] < -nlevel)
if not ilup.size and not ildown.size: if not ilup.size and not ildown.size:
print ("earllatepicker: Signal lower than noise level!") if stealth_mode is False:
print ("Skip this trace!") print ("earllatepicker: Signal lower than noise level!\n"
"Skip this trace!")
return LPick, EPick, PickError return LPick, EPick, PickError
il = min(np.min(ilup) if ilup.size else float('inf'), il = min(np.min(ilup) if ilup.size else float('inf'),
np.min(ildown) if ildown.size else float('inf')) np.min(ildown) if ildown.size else float('inf'))
@ -78,7 +80,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode=False):
# if EPick stays NaN the signal window size will be doubled # if EPick stays NaN the signal window size will be doubled
while np.isnan(EPick): while np.isnan(EPick):
if count > 0: if count > 0:
if stealthMode is False: if stealth_mode is False:
print("\nearllatepicker: Doubled signal window size %s time(s) " print("\nearllatepicker: Doubled signal window size %s time(s) "
"because of NaN for earliest pick." % count) "because of NaN for earliest pick." % count)
isigDoubleWinStart = pis[-1] + 1 isigDoubleWinStart = pis[-1] + 1
@ -87,7 +89,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode=False):
if (isigDoubleWinStart + len(pis)) < X[0].data.size: if (isigDoubleWinStart + len(pis)) < X[0].data.size:
pis = np.concatenate((pis, isignalDoubleWin)) pis = np.concatenate((pis, isignalDoubleWin))
else: else:
print("Could not double signal window. Index out of bounds.") if stealth_mode is False:
print("Could not double signal window. Index out of bounds.")
break break
count += 1 count += 1
# determine all zero crossings in signal window (demeaned) # determine all zero crossings in signal window (demeaned)

View File

@ -4,31 +4,151 @@
import os import os
import glob import glob
from obspy import UTCDateTime from obspy import UTCDateTime
import sys
def gse_time_header(lines):
'''
takes a path FILE to a GSE data file and returns the time header of the file
:param file: path to GSE data file
:type file: str
:return: time header from FILE
:rtype: str
>>> gse_time_header('test.gse')
"WID2 2005/10/09 20:17:25.000 ATH SHZ NSP CM6 9000 50.000000 0.10E+01 1.000 NSP -1.0 0.0"
'''
return lines[1]
def time_from_header(header): 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.
:param header: second line from .gse file
:type header: string
:return: a list of integers of form [year, month, day, hour, minute, second, microsecond]
"""
timeline = header.split(' ') timeline = header.split(' ')
time = timeline[1].split('/') + timeline[2].split(':') time = timeline[1].split('/') + timeline[2].split(':')
time = time[:-1] + time[-1].split('.') time = time[:-1] + time[-1].split('.')
time[-1] += '000'
return [int(t) for t in time] return [int(t) for t in time]
def check_time(time):
def check_time(datetime):
"""
Function takes in date and time as list and validates it's values by trying to make an UTCDateTime object from it
:param datetime: list of integers [year, month, day, hour, minute, second, microsecond]
:type datetime: list
:return: returns True if Values are in supposed range, returns False otherwise
>>> check_time([1999, 01, 01, 23, 59, 59, 999000])
True
>>> check_time([1999, 01, 01, 23, 59, 60, 999000])
False
>>> check_time([1999, 01, 01, 23, 59, 59, 1000000])
False
>>> check_time([1999, 01, 01, 23, 60, 59, 999000])
False
>>> check_time([1999, 01, 01, 23, 60, 59, 999000])
False
>>> check_time([1999, 01, 01, 24, 59, 59, 999000])
False
>>> check_time([1999, 01, 31, 23, 59, 59, 999000])
True
>>> check_time([1999, 02, 30, 23, 59, 59, 999000])
False
>>> check_time([1999, 02, 29, 23, 59, 59, 999000])
False
>>> check_time([2000, 02, 29, 23, 59, 59, 999000])
True
>>> check_time([2000, 13, 29, 23, 59, 59, 999000])
False
"""
try: try:
UTCDateTime(time) UTCDateTime(*datetime)
return True return True
except ValueError: except ValueError:
return False return False
def get_file_list(root_dir):
"""
Function uses a directorie to get all the *.gse files from it.
:param root_dir: a directorie leading to the .gse files
:type root_dir: string
:return: returns a list of filenames (without path to them)
"""
file_list = glob.glob1(root_dir, '*.gse')
return file_list
def checks_station_second(datetime, file):
"""
Function uses the given list to check if the parameter 'second' is set to 60 by mistake
and sets the time correctly if so. Can only correct time if no date change would be necessary.
:param datetime: [year, month, day, hour, minute, second, microsecond]
:return: returns the input with the correct value for second
"""
if datetime[5] == 60:
if datetime[4] == 59:
if datetime[3] == 23:
err_msg = 'Date should be next day. ' \
'File not changed: {0}'.format(file)
raise ValueError(err_msg)
else:
datetime[3] += 1
datetime[4] = 0
datetime[5] = 0
else:
datetime[4] += 1
datetime[5] = 0
return datetime
def make_time_line(line, datetime):
"""
Function takes in the original line from a .gse file and a list of date and
time values to make a new line with corrected date and time.
:param line: second line from .gse file.
:type line: string
:param datetime: list of integers [year, month, day, hour, minute, second, microsecond]
:type datetime: list
:return: returns a string to write it into a file.
"""
ins_form = '{0:02d}:{1:02d}:{2:02d}.{3:03d}'
insertion = ins_form.format(int(datetime[3]),
int(datetime[4]),
int(datetime[5]),
int(datetime[6] * 1e-3))
newline = line[:16] + insertion + line[28:]
return newline
def evt_head_check(root_dir, out_dir = None):
"""
A function to make sure that an arbitrary number of .gse files have correct values in their header.
:param root_dir: a directory leading to the .gse files.
:type root_dir: string
:param out_dir: a directory to store the new files somwhere els.
:return: returns nothing
"""
if not out_dir:
print('WARNING files are going to be overwritten!')
inp = str(raw_input('Continue? [y/N]'))
if not inp == 'y':
sys.exit()
filelist = get_file_list(root_dir)
nfiles = 0
for file in filelist:
infile = open(os.path.join(root_dir, file), 'r')
lines = infile.readlines()
infile.close()
datetime = time_from_header(lines[1])
if check_time(datetime):
continue
else:
nfiles += 1
datetime = checks_station_second(datetime, file)
print('writing ' + file)
# write File
lines[1] = make_time_line(lines[1], datetime)
if not out_dir:
out = open(os.path.join(root_dir, file), 'w')
out.writelines(lines)
out.close()
else:
out = open(os.path.join(out_dir, file), 'w')
out.writelines(lines)
out.close()
print(nfiles)
if __name__ == "__main__":
import doctest
doctest.testmod()

View File

@ -1,14 +1,14 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import os
import subprocess
import pwd
import re
import hashlib import hashlib
import numpy as np import numpy as np
import os
import pwd
import re
import subprocess
from obspy.core import UTCDateTime from obspy.core import UTCDateTime
import obspy.core.event as ope
def _pickle_method(m): def _pickle_method(m):
if m.im_self is None: if m.im_self is None:
@ -16,236 +16,20 @@ def _pickle_method(m):
else: else:
return getattr, (m.im_self, m.im_func.func_name) return getattr, (m.im_self, m.im_func.func_name)
def worker(func, input, cores = 'max', async = False):
def worker(func, input, cores='max', async=False):
return result
import multiprocessing import multiprocessing
if cores == 'max': if cores == 'max':
cores = multiprocessing.cpu_count() cores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(cores) pool = multiprocessing.Pool(cores)
if async == True: if async == True:
result = pool.map_async(func, input) result = pool.map_async(func, input)
else: else:
result = pool.map(func, input) result = pool.map(func, input)
pool.close() pool.close()
return result
def createAmplitude(pickID, amp, unit, category, cinfo):
'''
:param pickID:
:param amp:
:param unit:
:param category:
:param cinfo:
:return:
'''
amplitude = ope.Amplitude()
amplitude.creation_info = cinfo
amplitude.generic_amplitude = amp
amplitude.unit = ope.AmplitudeUnit(unit)
amplitude.type = ope.AmplitudeCategory(category)
amplitude.pick_id = pickID
return amplitude
def createArrival(pickresID, cinfo, phase, azimuth=None, dist=None):
'''
createArrival - function to create an Obspy Arrival
:param pickresID: Resource identifier of the created pick
:type pickresID: :class: `~obspy.core.event.ResourceIdentifier` object
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param phase: name of the arrivals seismic phase
:type phase: str
:param azimuth: azimuth between source and receiver
:type azimuth: float or int, optional
:param dist: distance between source and receiver
:type dist: float or int, optional
:return: An ObsPy :class: `~obspy.core.event.Arrival` object
'''
arrival = ope.Arrival()
arrival.creation_info = cinfo
arrival.pick_id = pickresID
arrival.phase = phase
if azimuth is not None:
arrival.azimuth = float(azimuth) if azimuth > -180 else azimuth + 360.
else:
arrival.azimuth = azimuth
arrival.distance = dist
return arrival
def createCreationInfo(agency_id=None, creation_time=None, author=None):
'''
:param agency_id:
:param creation_time:
:param author:
:return:
'''
if author is None:
author = getLogin()
if creation_time is None:
creation_time = UTCDateTime()
return ope.CreationInfo(agency_id=agency_id, author=author,
creation_time=creation_time)
def createEvent(origintime, cinfo, originloc=None, etype=None, resID=None,
authority_id=None):
'''
createEvent - funtion to create an ObsPy Event
:param origintime: the events origintime
:type origintime: :class: `~obspy.core.utcdatetime.UTCDateTime` object
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param originloc: tuple containing the location of the origin
(LAT, LON, DEP) affiliated with the event which is created
:type originloc: tuple, list
:param etype: Event type str object. converted via ObsPy to a valid event
type string.
:type etype: str
:param resID: Resource identifier of the created event
:type resID: :class: `~obspy.core.event.ResourceIdentifier` object, str
:param authority_id: name of the institution carrying out the processing
:type authority_id: str
:return: An ObsPy :class: `~obspy.core.event.Event` object
'''
etype = ope.EventType(etype)
if originloc is not None:
o = createOrigin(origintime, cinfo,
originloc[0], originloc[1], originloc[2])
else:
o = None
if etype is None:
etype = ope.EventType('earthquake') # defaults to 'earthquake'
if not resID:
resID = createResourceID(origintime, etype, authority_id)
elif isinstance(resID, str):
resID = createResourceID(origintime, etype, authority_id, resID)
elif not isinstance(resID, ope.ResourceIdentifier):
raise TypeError("unsupported type(resID) for resource identifier "
"generation: %s" % type(resID))
event = ope.Event(resource_id=resID)
event.creation_info = cinfo
event.event_type = etype
if o:
event.origins = [o]
return event
def createMagnitude(originID, cinfo):
'''
createMagnitude - function to create an ObsPy Magnitude object
:param originID:
:type originID:
:param cinfo:
:type cinfo:
:return:
'''
magnitude = ope.Magnitude()
magnitude.creation_info = cinfo
magnitude.origin_id = originID
return magnitude
def createOrigin(origintime, cinfo, latitude, longitude, depth):
'''
createOrigin - function to create an ObsPy Origin
:param origintime: the origins time of occurence
:type origintime: :class: `~obspy.core.utcdatetime.UTCDateTime` object
:param cinfo:
:type cinfo:
:param latitude: latitude in decimal degree of the origins location
:type latitude: float
:param longitude: longitude in decimal degree of the origins location
:type longitude: float
:param depth: hypocentral depth of the origin
:type depth: float
:return: An ObsPy :class: `~obspy.core.event.Origin` object
'''
assert isinstance(origintime, UTCDateTime), "origintime has to be " \
"a UTCDateTime object, but " \
"actually is of type " \
"'%s'" % type(origintime)
origin = ope.Origin()
origin.time = origintime
origin.creation_info = cinfo
origin.latitude = latitude
origin.longitude = longitude
origin.depth = depth
return origin
def createPick(origintime, picknum, picktime, eventnum, cinfo, phase, station,
wfseedstr, authority_id):
'''
createPick - function to create an ObsPy Pick
:param origintime:
:type origintime:
:param picknum: number of the created pick
:type picknum: int
:param picktime:
:type picktime:
:param eventnum: human-readable event identifier
:type eventnum: str
:param cinfo: An ObsPy :class: `~obspy.core.event.CreationInfo` object
holding information on the creation of the returned object
:type cinfo: :class: `~obspy.core.event.CreationInfo` object
:param phase: name of the arrivals seismic phase
:type phase: str
:param station: name of the station at which the seismic phase has been
picked
:type station: str
:param wfseedstr: A SEED formatted string of the form
network.station.location.channel in order to set a referenced waveform
:type wfseedstr: str, SEED formatted
:param authority_id: name of the institution carrying out the processing
:type authority_id: str
:return: An ObsPy :class: `~obspy.core.event.Pick` object
'''
pickID = eventnum + '_' + station.strip() + '/{0:03d}'.format(picknum)
pickresID = createResourceID(origintime, 'pick', authority_id, pickID)
pick = ope.Pick()
pick.resource_id = pickresID
pick.time = picktime
pick.creation_info = cinfo
pick.phase_hint = phase
pick.waveform_id = ope.ResourceIdentifier(id=wfseedstr, prefix='file:/')
return pick
def createResourceID(timetohash, restype, authority_id=None, hrstr=None):
'''
:param timetohash:
:type timetohash
:param restype: type of the resource, e.g. 'orig', 'earthquake' ...
:type restype: str
:param authority_id: name of the institution carrying out the processing
:type authority_id: str, optional
:param hrstr:
:type hrstr:
:return:
'''
assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \
"UTCDateTime object"
hid = getHash(timetohash)
if hrstr is None:
resID = ope.ResourceIdentifier(restype + '/' + hid[0:6])
else:
resID = ope.ResourceIdentifier(restype + '/' + hrstr + '_' + hid[0:6])
if authority_id is not None:
resID.convertIDToQuakeMLURI(authority_id=authority_id)
return resID
def demeanTrace(trace, window): def demeanTrace(trace, window):
@ -469,6 +253,7 @@ def runProgram(cmd, parameter=None):
output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd), output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd),
shell=True) shell=True)
if __name__ == "__main__": if __name__ == "__main__":
import doctest import doctest

View File

@ -19,15 +19,23 @@ if __name__ == '__main__':
) )
parser.add_argument( parser.add_argument(
'dbroot', type=str, help='specifies the root directory (in ' 'root', type=str, help='specifies the root directory'
'most cases PILOT database folder)'
) )
parser.add_argument( parser.add_argument(
'--output', '-o', type=str, help='path to the output directory', dest='output' 'db', type=str, help='specifies the database name'
) )
parser.add_argument( parser.add_argument(
'--parameterfile', '-p', type=str, help='full path to the parameterfile', dest='parfile' '--output', '-o', type=str, help='path to the output directory',
dest='output'
)
parser.add_argument(
'--parameterfile', '-p', type=str,
help='full path to the parameterfile', dest='parfile'
)
parser.add_argument(
'--verbosity', '-v', action='count', help='increase output verbosity',
default=0, dest='verbosity'
) )
args = parser.parse_args() args = parser.parse_args()
reassess_pilot_db(args.dbroot, args.output, args.parfile) reassess_pilot_db(args.root, args.db, args.output, args.parfile, args.verbosity)

View File

@ -19,8 +19,10 @@ if __name__ == '__main__':
) )
parser.add_argument( parser.add_argument(
'dbroot', type=str, help='specifies the root directory (in ' 'root', type=str, help='specifies the root directory'
'most cases PILOT database folder)' )
parser.add_argument(
'db', type=str, help='specifies the database name'
) )
parser.add_argument( parser.add_argument(
'id', type=str, help='PILOT event identifier' 'id', type=str, help='PILOT event identifier'
@ -33,4 +35,4 @@ if __name__ == '__main__':
) )
args = parser.parse_args() args = parser.parse_args()
reassess_pilot_event(args.dbroot, args.id, args.output, args.parfile) reassess_pilot_event(args.root, args.db, args.id, args.output, args.parfile)