diff --git a/QtPyLoT.py b/QtPyLoT.py index 11aa6719..6daaffa6 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -43,7 +43,7 @@ from obspy import UTCDateTime from pylot.core.io.data import Data from pylot.core.io.inputs import FilterOptions, AutoPickParameter from pylot.core.pick.autopick import autopickevent -from pylot.core.io.phases import picks_to_dict +from pylot.core.io.phases import picksdict_from_picks from pylot.core.loc.nll import locate as locateNll from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP,\ AUTOMATIC_DEFAULTS @@ -735,7 +735,7 @@ class MainWindow(QMainWindow): return rval def updatePicks(self, type='manual'): - picks = picks_to_dict(evt=self.getData().getEvtData()) + picks = picksdict_from_picks(evt=self.getData().getEvtData()) if type == 'manual': self.picks.update(picks) elif type == 'auto': diff --git a/pylot/core/io/data.py b/pylot/core/io/data.py index faa07cb8..4f869932 100644 --- a/pylot/core/io/data.py +++ b/pylot/core/io/data.py @@ -9,7 +9,7 @@ from obspy.core import read, Stream, UTCDateTime from obspy import read_events, read_inventory from obspy.core.event import Event, ResourceIdentifier, Pick, WaveformStreamID -from pylot.core.io.phases import readPILOTEvent, picks_from_dict +from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict from pylot.core.util.utils import fnConstructor, getGlobalTimes from pylot.core.util.errors import FormatError, OverwriteError @@ -411,19 +411,18 @@ class Data(object): :raise OverwriteError: raises an OverwriteError if the picks list is not empty. The GUI will then ask for a decision. """ - firstonset = None + + #firstonset = find_firstonset(picks) if self.getEvtData().picks: raise OverwriteError('Actual picks would be overwritten!') - picks, firstonset = picks_from_dict(picks) + picks = picks_from_picksdict(picks) self.getEvtData().picks = picks + # if 'smi:local' in self.getID() and firstonset: + # fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S') + # ID = ResourceIdentifier('event/' + fonset_str) + # ID.convertIDToQuakeMLURI(authority_id=authority_id) + # self.getEvtData().resource_id = ID - if 'smi:local' in self.getID() and firstonset: - fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S') - ID = ResourceIdentifier('event/' + fonset_str) - ID.convertIDToQuakeMLURI(authority_id=authority_id) - self.getEvtData().resource_id = ID - else: - print('No picks to apply!') def applyArrivals(arrivals): """ diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py index f3c6c96a..3924ddd3 100644 --- a/pylot/core/io/phases.py +++ b/pylot/core/io/phases.py @@ -8,9 +8,9 @@ import scipy.io as sio import obspy.core.event as ope from obspy.core import UTCDateTime +from pylot.core.pick.utils import select_for_phase from pylot.core.util.utils import getOwner, createPick, createArrival, \ - createEvent, createOrigin, createMagnitude - + createEvent, createOrigin, createMagnitude, getGlobalTimes def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): """ @@ -192,14 +192,14 @@ def picksdict_from_obs(fn): return picks -def picks_to_dict(evt): - ''' +def picksdict_from_picks(evt): + """ Takes an Event object and return the pick dictionary commonly used within PyLoT :param evt: Event object contain all available information :type evt: `~obspy.core.event.Event` :return: pick dictionary - ''' + """ picks = {} for pick in evt.picks: phase = {} @@ -229,12 +229,12 @@ def picks_to_dict(evt): picks[station] = onsets.copy() return picks -def picks_from_dict(picks): - firstonset = None +def picks_from_picksdict(picks): + 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): + if not isinstance(phase, dict) or len(phase) < 3: continue onset = phase['mpp'] epp = phase['epp'] @@ -262,52 +262,90 @@ def picks_from_dict(picks): else: pick.polarity = 'undecidable' except KeyError as e: - print('No polarity information found for %s' % phase) - if firstonset is None or firstonset > onset: - firstonset = onset + print(e.message, 'No polarity information found for %s' % phase) + picks_list.append(pick) + return picks_list -def reassess_pilot_event(root_dir, event_id): +def reassess_pilot_db(root_dir, out_dir=None, fn_param=None): + import glob + + evt_list = glob.glob1(root_dir,'e????.???.??') + + for evt in evt_list: + reassess_pilot_event(root_dir, evt, out_dir, fn_param) + + +def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbosity=0): from obspy import read - from pylot.core.util.defaults import AUTOMATIC_DEFAULTS + from pylot.core.io.inputs import AutoPickParameter from pylot.core.pick.utils import earllatepicker - default = AutoPickParameter(AUTOMATIC_DEFAULTS) + if fn_param is None: + import pylot.core.util.defaults as defaults + fn_param = defaults.AUTOMATIC_DEFAULTS + + default = AutoPickParameter(fn_param) search_base = os.path.join(root_dir, event_id) phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat')) - picks_dict = picks_from_pilot(phases_file) + if not phases_file: + return + print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0])) + picks_dict = picks_from_pilot(phases_file[0]) + print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict)) for station in picks_dict.keys(): fn_pattern = os.path.join(search_base, '{0}*'.format(station)) try: - st = read(fn_pattern) + st = read(fn_pattern, format='GSE2') except TypeError as e: print(e.message) - st = read(fn_pattern, format='GSE2') - if not st: - raise RuntimeError('no waveform data found for station {station}'.format(station=station)) + st = read(fn_pattern) for phase in picks_dict[station].keys(): try: mpp = picks_dict[station][phase]['mpp'] except KeyError as e: print(e.message, station) continue - epp, lpp, spe = earllatepicker(st, + sel_st = select_for_phase(st, phase) + if not sel_st: + raise warnings.formatwarning( + 'no waveform data found for station {station}'.format( + station=station), category=RuntimeWarning) + print(sel_st) + stime, etime = getGlobalTimes(sel_st) + rel_pick = mpp - stime + epp, lpp, spe = earllatepicker(sel_st, default.get('nfac{0}'.format(phase)), default.get('tsnrz' if phase == 'P' else 'tsnrh'), - mpp + Pick1=rel_pick, + iplot=None, ) + if epp is None or lpp is None: + continue + epp = stime + epp + lpp = stime + lpp + min_diff = 3 * st[0].stats.delta + if lpp - mpp < min_diff: + lpp = mpp + min_diff + if mpp - epp < min_diff: + epp = mpp - min_diff picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe) # create Event object for export evt = ope.Event(resource_id=event_id) - evt.picks = picks_from_dict(picks_dict) + evt.picks = picks_from_picksdict(picks_dict) # write phase information to file - evt.write('{0}.xml'.format(event_id), format='QUAKEML') + if not out_dir: + fnout_prefix = os.path.join(root_dir, event_id, '{0}.'.format(event_id)) + else: + fnout_prefix = os.path.join(out_dir, '{0}.'.format(event_id)) + evt.write(fnout_prefix + 'xml', format='QUAKEML') + #evt.write(fnout_prefix + 'cnv', format='VELEST') def writephases(arrivals, fformat, filename): - ''' + """ Function of methods to write phases to the following standard file formats used for locating earthquakes: @@ -325,7 +363,7 @@ def writephases(arrivals, fformat, filename): :param: filename, full path and name of phase file :type: string - ''' + """ if fformat == 'NLLoc': print ("Writing phases to %s for NLLoc" % filename) @@ -387,7 +425,6 @@ def writephases(arrivals, fformat, filename): sweight)) fid.close() - elif fformat == 'HYPO71': print ("Writing phases to %s for HYPO71" % filename) fid = open("%s" % filename, 'w') diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py index d4927ead..1a828c2a 100644 --- a/pylot/core/pick/compare.py +++ b/pylot/core/pick/compare.py @@ -9,7 +9,7 @@ import matplotlib.pyplot as plt from obspy import read_events -from pylot.core.io.phases import picks_to_dict +from pylot.core.io.phases import picksdict_from_picks from pylot.core.util.pdf import ProbabilityDensityFunction from pylot.core.util.version import get_git_version as _getVersionString @@ -86,8 +86,8 @@ class Comparison(object): """ compare_pdfs = dict() - pdf_a = self.get(self.names[0]).pdf_data(type) - pdf_b = self.get(self.names[1]).pdf_data(type) + pdf_a = self.get(self.names[0]).generate_pdf_data(type) + pdf_b = self.get(self.names[1]).generate_pdf_data(type) for station, phases in pdf_a.items(): if station in pdf_b.keys(): @@ -206,6 +206,7 @@ class PDFDictionary(object): def __init__(self, data): self._pickdata = data + self._pdfdata = self.generate_pdf_data() def __nonzero__(self): if len(self.pick_data) < 1: @@ -213,6 +214,14 @@ class PDFDictionary(object): else: return True + @property + def pdf_data(self): + return self._pdfdata + + @pdf_data.setter + def pdf_data(self, data): + self._pdfdata = data + @property def pick_data(self): return self._pickdata @@ -221,15 +230,23 @@ class PDFDictionary(object): def pick_data(self, data): self._pickdata = data + @property + def stations(self): + return self.pick_data.keys() + + @property + def nstations(self): + return len(self.stations) + @classmethod def from_quakeml(self, fn): cat = read_events(fn) if len(cat) > 1: raise NotImplementedError('reading more than one event at the same ' 'time is not implemented yet! Sorry!') - return PDFDictionary(picks_to_dict(cat[0])) + return PDFDictionary(picksdict_from_picks(cat[0])) - def pdf_data(self, type='exp'): + def generate_pdf_data(self, type='exp'): """ Returns probabiliy density function dictionary containing the representation of the actual pick_data. @@ -251,3 +268,57 @@ class PDFDictionary(object): return pdf_picks + def plot(self, stations=None): + ''' + plots the all probability density function for either desired STATIONS + or all available date + :param stations: list of stations to be plotted + :type stations: list + :return: matplotlib figure object containing the plot + ''' + assert stations is not None or not isinstance(stations, list), \ + 'parameter stations should be a list not {0}'.format(type(stations)) + if not stations: + nstations = self.nstations + stations = self.stations + else: + nstations = len(stations) + + istations = range(nstations) + fig, axarr = plt.subplots(nstations, 2, sharex='col', sharey='row') + hide_labels = True + + for n in istations: + station = stations[n] + pdfs = self.pdf_data[station] + for l, phase in enumerate(pdfs.keys()): + try: + axarr[n, l].plot(pdfs[phase].axis, pdfs[phase].data) + if n is 0: + axarr[n, l].set_title(phase) + if l is 0: + axann = axarr[n, l].annotate(station, xy=(.05, .5), + xycoords='axes fraction') + bbox_props = dict(boxstyle='round', facecolor='lightgrey', + alpha=.7) + axann.set_bbox(bbox_props) + if n == int(np.median(istations)) and l is 0: + label = 'probability density (qualitative)' + axarr[n, l].set_ylabel(label) + except IndexError as e: + print('trying aligned plotting\n{0}'.format(e)) + hide_labels = False + axarr[l].plot(pdfs[phase].axis, pdfs[phase].data) + axarr[l].set_title(phase) + if l is 0: + axann = axarr[l].annotate(station, xy=(.05, .5), + xycoords='axes fraction') + bbox_props = dict(boxstyle='round', facecolor='lightgrey', + alpha=.7) + axann.set_bbox(bbox_props) + if hide_labels: + plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False) + plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False) + plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False) + + return fig diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index b8cf358d..e603cfa4 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -464,6 +464,37 @@ def getResolutionWindow(snr): return time_resolution / 2 +def select_for_phase(st, phase): + ''' + takes a STream object and a phase name and returns that particular component + which presumably shows the chosen PHASE best + + :param st: stream object containing one or more component[s] + :type st: `~obspy.core.stream.Stream` + :param phase: label of the phase for which the stream selection is carried + out; 'P' or 'S' + :type phase: str + :return: + ''' + from pylot.core.util.defaults import COMPNAME_MAP + + sel_st = Stream() + if phase.upper() == 'P': + comp = 'Z' + alter_comp = COMPNAME_MAP[comp] + sel_st += st.select(component=comp) + sel_st += st.select(component=alter_comp) + elif phase.upper() == 'S': + comps = 'NE' + for comp in comps: + alter_comp = COMPNAME_MAP[comp] + sel_st += st.select(component=comp) + sel_st += st.select(component=alter_comp) + else: + raise TypeError('Unknown phase label: {0}'.format(phase)) + return sel_st + + def wadaticheck(pickdic, dttolerance, iplot): ''' Function to calculate Wadati-diagram from given P and S onsets in order @@ -938,10 +969,6 @@ def checkZ4S(X, pick, zfac, checkwin, iplot): return returnflag -def reassess_pilot_event(): - pass - - if __name__ == '__main__': import doctest diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py new file mode 100644 index 00000000..f91e2ba2 --- /dev/null +++ b/pylot/core/util/dataprocessing.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +import glob +from obspy import UTCDateTime + +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): + timeline = header.split(' ') + time = timeline[1].split('/') + timeline[2].split(':') + time = time[:-1] + time[-1].split('.') + time[-1] += '000' + return [int(t) for t in time] + +def check_time(time): + try: + UTCDateTime(time) + return True + except ValueError: + return False \ No newline at end of file diff --git a/pylot/core/util/pdf.py b/pylot/core/util/pdf.py index 776192aa..1a233aed 100644 --- a/pylot/core/util/pdf.py +++ b/pylot/core/util/pdf.py @@ -256,9 +256,10 @@ class ProbabilityDensityFunction(object): ''' rval = 0 - for n, x in enumerate(self.axis): + axis = self.axis - self.x0 + for n, x in enumerate(axis): rval += x * self.data[n] - return rval * self.incr + return rval * self.incr + self.x0 def standard_deviation(self): mu = self.expectation() diff --git a/scripts/pylot-noisewindow.py b/scripts/pylot-noisewindow.py old mode 100644 new mode 100755 diff --git a/scripts/pylot-reasses-pilot-db.py b/scripts/pylot-reasses-pilot-db.py new file mode 100755 index 00000000..5bb4441d --- /dev/null +++ b/scripts/pylot-reasses-pilot-db.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import argparse + +from pylot.core.util.version import get_git_version as _getVersionString +from pylot.core.io.phases import reassess_pilot_db + +__version__ = _getVersionString() +__author__ = 'S. Wehling-Benatelli' + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description='reassess old PILOT event data base in terms of consistent ' + 'automatic uncertainty estimation', + epilog='Script written by {author} belonging to PyLoT version' + ' {version}\n'.format(author=__author__, + version=__version__) + ) + + parser.add_argument( + 'dbroot', type=str, help='specifies the root directory (in ' + 'most cases PILOT database folder)' + ) + parser.add_argument( + '--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' + ) + + args = parser.parse_args() + reassess_pilot_db(args.dbroot, args.output, args.parfile) diff --git a/scripts/pylot-reasses-pilot-event.py b/scripts/pylot-reasses-pilot-event.py old mode 100644 new mode 100755 index 074f5218..1852e680 --- a/scripts/pylot-reasses-pilot-event.py +++ b/scripts/pylot-reasses-pilot-event.py @@ -3,23 +3,34 @@ import argparse -from pylot.core.pick.utils import reassess_pilot_event from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.io.phases import reassess_pilot_event __version__ = _getVersionString() -__author__ = 'sebastianw' +__author__ = 'S. Wehling-Benatelli' if __name__ == '__main__': - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser( + description='reassess old PILOT event data in terms of consistent ' + 'automatic uncertainty estimation', + epilog='Script written by {author} belonging to PyLoT version' + ' {version}\n'.format(author=__author__, + version=__version__) + ) parser.add_argument( - '--directory', '-d', type=str, help='specifies the root directory (in ' + 'dbroot', type=str, help='specifies the root directory (in ' 'most cases PILOT database folder)' ) parser.add_argument( - '--eventid', '-i', type=str, help='PILOT event identifier' + 'id', type=str, help='PILOT event identifier' + ) + parser.add_argument( + '--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' ) args = parser.parse_args() - reassess_pilot_event(args.dir, args.id) + reassess_pilot_event(args.dbroot, args.id, args.output, args.parfile) diff --git a/scripts/pylot-signalwindow.py b/scripts/pylot-signalwindow.py old mode 100644 new mode 100755 diff --git a/scripts/pylot-snr.py b/scripts/pylot-snr.py old mode 100644 new mode 100755