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

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

Browse Source
This commit is contained in:
Marcel Paffrath
2016-05-31 11:15:38 +02:00
parent 0520df5963 51aa949b56
commit e341994825
7 changed files with 239 additions and 166 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -1,17 +1,17 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import glob
from obspy.io.xseed import Parser
from obspy.core import read, Stream, UTCDateTime
import os
from obspy import read_events, read_inventory
from obspy.core.event import Event, ResourceIdentifier, Pick, WaveformStreamID
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
from pylot.core.util.utils import fnConstructor, getGlobalTimes
from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \
picksdict_from_pilot
from pylot.core.util.errors import FormatError, OverwriteError
from pylot.core.util.utils import fnConstructor, getGlobalTimes
class Data(object):
@@ -36,17 +36,36 @@ class Data(object):
self.wfdata = Stream()
self._new = False
if isinstance(evtdata, Event):
self.evtdata = evtdata
pass
elif isinstance(evtdata, dict):
evt = readPILOTEvent(**evtdata)
self.evtdata = evt
elif evtdata:
cat = read_events(evtdata)
self.evtdata = cat[0]
evtdata = evt
elif isinstance(evtdata, str):
try:
cat = read_events(evtdata)
if len(cat) is not 1:
raise ValueError('ambiguous event information for file: '
'{file}'.format(file=evtdata))
evtdata = cat[0]
except TypeError as e:
if 'Unknown format for file' in e.message:
if 'PHASES' in evtdata:
picks = picksdict_from_pilot(evtdata)
evtdata = Event()
evtdata.picks = picks_from_picksdict(picks)
elif 'LOC' in evtdata:
raise NotImplementedError('PILOT location information '
'read support not yet '
'implemeted.')
else:
raise e
else:
raise e
else: # create an empty Event object
self.setNew()
self.evtdata = Event()
self.getEvtData().picks = []
evtdata = Event()
evtdata.picks = []
self.evtdata = evtdata
self.wforiginal = None
self.cuttimes = None
self.dirty = False

10
pylot/core/io/location.py
View File

@@ -68,8 +68,8 @@ def create_creation_info(agency_id=None, creation_time=None, author=None):
creation_time=creation_time)
def create_event(origintime, cinfo, originloc=None, etype=None, resID=None,
authority_id=None):
def create_event(origintime, cinfo, originloc=None, etype='earthquake',
resID=None, authority_id=None):
'''
create_event - funtion to create an ObsPy Event
@@ -90,14 +90,12 @@ def create_event(origintime, cinfo, originloc=None, etype=None, resID=None,
: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):
@@ -216,7 +214,7 @@ def create_resourceID(timetohash, restype, authority_id=None, hrstr=None):
if hrstr is None:
resID = ope.ResourceIdentifier(restype + '/' + hid[0:6])
else:
resID = ope.ResourceIdentifier(restype + '/' + hrstr + '_' + hid[0:6])
resID = ope.ResourceIdentifier(restype + '/' + hrstr)
if authority_id is not None:
resID.convertIDToQuakeMLURI(authority_id=authority_id)
return resID

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

@@ -8,13 +8,14 @@ import scipy.io as sio
import warnings
from obspy.core import UTCDateTime
from pylot.core.io.inputs import AutoPickParameter
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.util.utils import getOwner, getGlobalTimes
from pylot.core.util.utils import getOwner, getGlobalTimes, four_digits
def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
"""
readPILOTEvent - function
@@ -30,7 +31,6 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
:param locfn: filename of the old PILOT Matlab LOC file
:return event: event object containing event and phase information
"""
sdir = os.path.split(phasfn)[0]
if phasfn is not None and os.path.isfile(phasfn):
phases = sio.loadmat(phasfn)
phasctime = UTCDateTime(os.path.getmtime(phasfn))
@@ -53,96 +53,55 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
loccinfo = ope.CreationInfo(agency_id=authority_id,
author=locauthor,
creation_time=locctime)
np = 0
try:
eventNum = loc['ID'][0]
# retrieve eventID for the actual database
idsplit = eventNum.split('.')
eventNum = str(loc['ID'][0])
# retrieve date information
julday = int(idsplit[1])
year = int(idsplit[2])
hour = int(loc['hh'])
minute = int(loc['mm'])
second = int(loc['ss'])
# retrieve eventID for the actual database
idsplit = eventNum.split('.')
if year + 2000 < UTCDateTime.utcnow().year:
year += 2000
else:
year += 1900
# retrieve date information
julday = int(idsplit[1])
year = int(idsplit[2])
hour = int(loc['hh'])
minute = int(loc['mm'])
second = int(loc['ss'])
eventDate = UTCDateTime(year=year, julday=julday, hour=hour,
minute=minute, second=second)
year = four_digits(year)
stations = [stat for stat in phases['stat'][0:-1:3]]
eventDate = UTCDateTime(year=year, julday=julday, hour=hour,
minute=minute, second=second)
event = create_event(eventDate, loccinfo, etype='earthquake', resID=eventNum,
authority_id=authority_id)
stations = [stat for stat in phases['stat'][0:-1:3]]
lat = float(loc['LAT'])
lon = float(loc['LON'])
dep = float(loc['DEP'])
lat = float(loc['LAT'])
lon = float(loc['LON'])
dep = float(loc['DEP'])
origin = create_origin(eventDate, loccinfo, lat, lon, dep)
for n, pick in enumerate(phases['Ptime']):
if pick[0] > 0:
kwargs = {'year': int(pick[0]),
'month': int(pick[1]),
'day': int(pick[2]),
'hour': int(pick[3]),
'minute': int(pick[4]),
'second': int(str(pick[5]).split('.')[0]),
'microsecond': int(str(pick[5]).split('.')[1][0:6])}
spick = phases['Stime'][n]
if spick[0] > 0:
skwargs = {'year': int(spick[0]),
'month': int(spick[1]),
'day': int(spick[2]),
'hour': int(spick[3]),
'minute': int(spick[4]),
'second': int(str(spick[5]).split('.')[0]),
'microsecond': int(str(spick[5]).split('.')[1][0:6])}
spicktime = UTCDateTime(**skwargs)
else:
spicktime = None
ppicktime = UTCDateTime(**kwargs)
event = create_event(eventDate, loccinfo, originloc=(lat, lon, dep),
etype='earthquake', resID=eventNum,
authority_id=authority_id)
for picktime, phase in [(ppicktime, 'P'), (spicktime, 'S')]:
if picktime is not None:
if phase == 'P':
wffn = os.path.join(sdir, '{0}*{1}*'.format(
stations[n].strip(), 'z'))
else:
wffn = os.path.join(sdir, '{0}*{1}*'.format(
stations[n].strip(), '[ne]'))
print(wffn)
pick = create_pick(eventDate, np, picktime, eventNum, pickcinfo,
phase, stations[n], wffn, authority_id)
event.picks.append(pick)
pickID = pick.get('id')
arrival = create_arrival(pickID, pickcinfo, phase)
origin.arrivals.append(arrival)
np += 1
picks = picksdict_from_pilot(phasfn)
event.picks = picks_from_picksdict(picks, creation_info=pickcinfo)
if event.origins:
origin = event.origins[0]
magnitude = create_magnitude(origin.get('id'), loccinfo)
magnitude.mag = float(loc['Mnet'])
magnitude.magnitude_type = 'Ml'
event.picks.append(pick)
event.origins.append(origin)
event.magnitudes.append(magnitude)
return event
except AttributeError as e:
raise AttributeError('{0} - Matlab LOC files {1} and {2} contains \
insufficient data!'.format(e, phasfn, locfn))
return event
def picks_from_pilot(fn):
def picksdict_from_pilot(fn):
from pylot.core.util.defaults import TIMEERROR_DEFAULTS
picks = dict()
phases_pilot = sio.loadmat(fn)
stations = stations_from_pilot(phases_pilot['stat'])
params = AutoPickParameter(TIMEERROR_DEFAULTS)
timeerrors = dict(P=params.get('timeerrorsP'),
S=params.get('timeerrorsS'))
for n, station in enumerate(stations):
phases = dict()
for onset_name in 'PS':
@@ -151,7 +110,14 @@ def picks_from_pilot(fn):
if not pick[0]:
continue
pick = convert_pilot_times(pick)
phases[onset_name] = dict(mpp=pick)
uncertainty_label = '{0}weight'.format(onset_name.lower())
ierror = phases_pilot[uncertainty_label][0, n]
try:
spe = timeerrors[onset_name][ierror]
except IndexError as e:
print(e.message + '\ntake two times the largest default error value')
spe = timeerrors[onset_name][-1] * 2
phases[onset_name] = dict(mpp=pick, spe=spe)
picks[station] = phases
return picks
@@ -231,27 +197,31 @@ def picksdict_from_picks(evt):
picks[station] = onsets.copy()
return picks
def picks_from_picksdict(picks):
def picks_from_picksdict(picks, creation_info=None):
picks_list = list()
for station, onsets in picks.items():
#print('Reading picks on station %s' % station)
for label, phase in onsets.items():
if not isinstance(phase, dict) or len(phase) < 3:
if not isinstance(phase, dict):
continue
onset = phase['mpp']
epp = phase['epp']
lpp = phase['lpp']
pick = ope.Pick()
if creation_info:
pick.creation_info = creation_info
pick.time = onset
error = phase['spe']
pick.time_errors.uncertainty = error
try:
epp = phase['epp']
lpp = phase['lpp']
pick.time_errors.lower_uncertainty = onset - epp
pick.time_errors.upper_uncertainty = lpp - onset
except KeyError as e:
warnings.warn(e.message, RuntimeWarning)
try:
picker = phase['picker']
except KeyError as e:
warnings.warn(str(e), Warning)
warnings.warn(e.message, RuntimeWarning)
picker = 'Unknown'
pick = ope.Pick()
pick.time = onset
pick.time_errors.lower_uncertainty = onset - epp
pick.time_errors.upper_uncertainty = lpp - onset
pick.time_errors.uncertainty = error
pick.phase_hint = label
pick.method_id = ope.ResourceIdentifier(id=picker)
pick.waveform_id = ope.WaveformStreamID(station_code=station)
@@ -303,7 +273,7 @@ def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None
return
if verbosity > 1:
print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0]))
picks_dict = picks_from_pilot(phases_file[0])
picks_dict = picksdict_from_pilot(phases_file[0])
if verbosity > 0:
print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict))
datacheck = list()

4
pylot/core/util/defaults.py
View File

@@ -47,6 +47,10 @@ AUTOMATIC_DEFAULTS = os.path.join(os.path.expanduser('~'),
'.pylot',
'autoPyLoT.in')
TIMEERROR_DEFAULTS = os.path.join(os.path.expanduser('~'),
'.pylot',
'PILOT_TimeErrors.in')
OUTPUTFORMATS = {'.xml': 'QUAKEML',
'.cnv': 'CNV',
'.obs': 'NLLOC_OBS'}

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

@@ -22,8 +22,8 @@ def worker(func, input, cores='max', async=False):
if cores == 'max':
cores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(cores)
pool = multiprocessing.Pool(cores)
if async == True:
result = pool.map_async(func, input)
else:
@@ -93,6 +93,14 @@ def fnConstructor(s):
return fn
def four_digits(year):
if year + 2000 < UTCDateTime.utcnow().year:
year += 2000
else:
year += 1900
return year
def getGlobalTimes(stream):
'''