From 46a20a10e687306848051f74c5753bbda397fb3d Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 9 Feb 2015 13:24:55 +0100 Subject: [PATCH] new functions added for event creation purposes: getHash - returns a hash string from an UTCDateTime object createResourceID - returns a valid PyLoT resourceID for arbitrary types of event data createOrigin - returns an ObsPy Origin object (work in progress) createEvent - returns an ObsPy Event object (work in progress) createPick - returns an ObsPy Pick object (work in progress) createArrival - returns an ObsPy Arrival object (work in progress) createMagnitude - returns an ObsPy Magnitude object (work in progress) createAmplitude - returns an ObsPy Amplitude object (work in progress) testing should be carried out as a next step --- pylot/core/read/io.py | 131 +++++++++-------------- pylot/core/util/__init__.py | 17 +-- pylot/core/util/utils.py | 205 +++++++++++++++++++++++++++++++++--- 3 files changed, 247 insertions(+), 106 deletions(-) diff --git a/pylot/core/read/io.py b/pylot/core/read/io.py index ff6ee866..3147ebe2 100644 --- a/pylot/core/read/io.py +++ b/pylot/core/read/io.py @@ -2,11 +2,14 @@ # -*- coding: utf-8 -*- import os -import glob + import scipy.io as sio import obspy.core.event as ope from obspy.core import UTCDateTime -from pylot.core.util import getOwner, createPick, createArrival, createEvent + +from pylot.core.util import getOwner, createPick, createArrival, createEvent, \ + createOrigin, createMagnitude + def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): """ @@ -61,7 +64,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): minute = int(loc['mm']) second = int(loc['ss']) - if UTCDateTime(year=year+2000) < UTCDateTime.utcnow(): + if UTCDateTime(year=year + 2000) < UTCDateTime.utcnow(): year += 2000 else: year += 1900 @@ -71,30 +74,31 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): stations = [stat for stat in phases['stat'][0:-1:3]] - eventid = 'loc/' + eventNum - evresID = ope.ResourceIdentifier(id=eventid) - evresID.convertIDToQuakeMLURI(authority_id=authority_id) - event = ope.Event(resource_id=evresID) - event.creation_info = loccinfo - etype = ope.EventType('earthquake') - event.event_type = etype + event = createEvent(eventDate, loccinfo, 'earthquake', eventNum, + authority_id) + + lat = float(loc['LAT']) + lon = float(loc['LON']) + dep = float(loc['DEP']) + + origin = createOrigin(eventDate, loccinfo, lat, lon, dep, eventNum) for n, pick in enumerate(phases['Ptime']): - 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])} + 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])} + 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 @@ -102,69 +106,34 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): for picktime, phase in [(ppicktime, 'P'), (spicktime, 'S')]: if phase == 'P': - wffn = os.path.join([sdir, '{0}*{1}*'.format(stations[n], 'z')]) + wffn = os.path.join([sdir, '{0}*{1}*'.format(stations[n], + 'z')]) else: - wffn = os.path.join([sdir, '{0}*{1}*'.format(stations[n], '[ne]')]) - pick, arrival, np = createArrival(np, picktime, eventNum, - stations[n], pickcinfo, phase, - wffn, authority_id) + wffn = os.path.join([sdir, '{0}*{1}*'.format(stations[n], + '[ne]')]) + pick = createPick(eventDate, np, picktime, eventNum, pickcinfo, + phase, + stat[n], wffn, authority_id) event.picks.append(pick) + pickID = pick.get('resource_id') + arrival = createArrival(pickID, eventNum, pickcinfo, phase, + stat[n], authority_id) + origin.arrivals.append(arrival) + np += 1 - origID = 'orig/' + genID - origresID = ResourceIdentifier(id=origID) - origresID.convertIDToQuakeMLURI(authority_id='BUG') - origin = Origin() - origin.resource_id = origresID - otime = self.location[eventid]['Location']['Origin time'] - origin.time = UTCDateTime(otime) - origin.creation_info = self.cinfo - HW = self.location[eventid]['Location']['Hochwert'] - RW = self.location[eventid]['Location']['Rechtswert'] - LAT, LON = coordTrans(HW, RW, tosys=GPS) - origin.latitude = LAT - origin.longitude = LON - origin.depth = 1000. - origin.depth_type = OriginDepthType('operator assigned') - origin.evaluation_mode = 'automatic' - origin.arrivals.append(arrival) + magnitude = createMagnitude(origin.get('resource_id'), eventDate, + loccinfo, + authority_id) + magnitude.mag = float(loc['Mnet']) + magnitude.magnitude_type = 'Ml' - amplID = 'corrampl/' + genID - amplresID = ResourceIdentifier(id=amplID) - amplresID.convertIDToQuakeMLURI(authority_id='BUG') - amplitude = Amplitude() - amplitude.resource_id = amplresID - amplitude.creation_info = self.cinfo - amp = self.data[eventid][phase]['Amplitude']*1e-9 - amplitude.generic_amplitude = amp - amplitude.unit = AmplitudeUnit('m/s') - amplitude.magnitude_hint = 'Ml' - amplitude.type = AmplitudeCategory('point') - amplitude.pick_id = pickresID + event.picks.append(pick) + event.origins.append(origin) + event.magnitudes.append(magnitude) - magnID = 'corrmag/' + genID - magnresID = ResourceIdentifier(id=magnID) - magnresID.convertIDToQuakeMLURI(authority_id='BUG') - magnitude = Magnitude() - magnitude.resource_id = magnresID - magnitude.creation_info = self.cinfo - magnitude.origin_id = origresID - mag = self.location[eventid]['Ml'] - magnitude.mag = mag - magnitude.magnitude_type = 'Ml' - - event = Event(resource_id=evresID) - event.creation_info = self.cinfo - etype = EventType('earthquake') - event.event_type = etype - edesc = EventDescription(text='Prosper Haniel (induced)') - event.event_descriptions.append(edesc) - event.picks.append(pick) - event.origins.append(origin) - event.magnitudes.append(magnitude) - event.amplitudes.append(amplitude) except AttributeError, e: - raise AttributeError('{0} - Matlab LOC file contain \ - insufficient data!'.format(e)) + raise AttributeError('{0} - Matlab LOC files {1} and {2} contains \ + insufficient data!'.format(e, phasfn, locfn)) diff --git a/pylot/core/util/__init__.py b/pylot/core/util/__init__.py index 4a856b7a..88c0ec76 100755 --- a/pylot/core/util/__init__.py +++ b/pylot/core/util/__init__.py @@ -1,18 +1,11 @@ from pylot.core.util.connection import checkurl from pylot.core.util.defaults import FILTERDEFAULTS -from pylot.core.util.errors import OptionsError -from pylot.core.util.errors import FormatError +from pylot.core.util.errors import OptionsError, FormatError from pylot.core.util.layouts import layoutStationButtons -from pylot.core.util.utils import fnConstructor -from pylot.core.util.utils import createEvent -from pylot.core.util.utils import getOwner -from pylot.core.util.utils import createArrival -from pylot.core.util.widgets import PickDlg -from pylot.core.util.widgets import HelpForm -from pylot.core.util.widgets import FilterOptionsDialog -from pylot.core.util.widgets import PropertiesDlg -from pylot.core.util.widgets import NewEventDlg -from pylot.core.util.widgets import MPLWidget +from pylot.core.util.utils import fnConstructor, createArrival, createEvent,\ + createPick, createOrigin, createMagnitude, getOwner, getHash +from pylot.core.util.widgets import PickDlg, HelpForm, FilterOptionsDialog,\ + PropertiesDlg, NewEventDlg, MPLWidget from pylot.core.util.version import get_git_version as _getVersionString diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 1ec04f60..583a0cbb 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -5,10 +5,12 @@ import os import pwd import re +import hashlib +from obspy.core import UTCDateTime import obspy.core.event as ope -def fnConstructor(s): +def fnConstructor(s): s = s.split('/')[-1] badchars = re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$') @@ -20,21 +22,161 @@ def fnConstructor(s): fn = '_' + fn return fn -def createEvent(origintime, latitude, longitude, depth, **kwargs): - evt = ope.Event() -def createArrival(picknum, picktime, eventnum, station, cinfo, phase, wfname, - authority_id): - pickID = 'pick/' + eventnum + '/' + station + '/{0:3d}'.format(picknum) - pickresID = ope.ResourceIdentifier(id=pickID) - pickresID.convertIDToQuakeMLURI(authority_id=authority_id) +def getHash(time): + ''' + + :param time: time object for which a hash should be calculated + :type time: :class: `~obspy.core.utcdatetime.UTCDateTime` object + :return: str + ''' + hg = hashlib.sha1() + hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f')) + return hg.hexdigest() + + +def createResourceID(timetohash, restype, authority_id=None, hrstr=None): + ''' + + :param 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 + :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 createOrigin(origintime, cinfo, latitude, longitude, depth, resID=None, + authority_id=None): + ''' + createOrigin - function to create an ObsPy Origin + :param origintime: the origins time of occurence + :type origintime: :class: `~obspy.core.utcdatetime.UTCDateTime` object + :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 + ''' + if resID is None: + resID = createResourceID(origintime, 'orig', authority_id=authority_id) + elif isinstance(resID, str): + resID = createResourceID(origintime, 'orig', authority_id=authority_id, + hrstr=resID) + + origin = ope.Origin() + origin.resource_id = resID + origin.time = UTCDateTime(origintime) + origin.creation_info = cinfo + origin.latitude = latitude + origin.longitude = longitude + origin.depth = depth + return origin + + +def createEvent(origintime, cinfo, etype, 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 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 + :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 etype is None: + etype = ope.EventType('earthquake') # defaults to 'earthquake' + if resID is None: + resID = createResourceID(origintime, etype, authority_id) + elif isinstance(resID, str): + resID = createResourceID(origintime, etype, authority_id, resID) + event = ope.Event(resource_id=resID) + event.creation_info = cinfo + event.event_type = etype + return event + + +def createPick(origintime, picknum, picktime, eventnum, cinfo, phase, station, + wfseedstr, + authority_id): + ''' + createPick - function to create an ObsPy Pick + + :param picknum: number of the created pick + :type picknum: int + :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 + '/{0:3d}'.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=wfname, prefix='file:/') + pick.waveform_id = ope.ResourceIdentifier(id=wfseedstr, prefix='file:/') + return pick + +def createArrival(pickresID, eventnum, cinfo, phase, station, authority_id, + 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 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 authority_id: name of the institution carrying out the processing + :type authority_id: 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 + ''' arriID = 'arrival/' + eventnum + '/' + station + '/{0}'.format(phase) arriresID = ope.ResourceIdentifier(id=arriID) arriresID.convertIDToQuakeMLURI(authority_id=authority_id) @@ -42,10 +184,47 @@ def createArrival(picknum, picktime, eventnum, station, cinfo, phase, wfname, arrival.resource_id = arriresID arrival.creation_info = cinfo arrival.pick_id = pickresID - arrival.phase = pick.phase_hint - azi = self.location[eventid]['Backazimuth'] - 180 - arrival.azimuth = azi if azi > -180 else azi + 360 - arrival.distance = self.location[eventid]['Distance']['deg'] + arrival.phase = phase + if azimuth is not None: + arrival.azimuth = float(azimuth) if azimuth > -180 else azimuth + 360 + else: + arrival.azimuth = azimuth + arrival.distance = None + return arrival + + +def createMagnitude(originID, origintime, cinfo, authority_id=None): + ''' + createMagnitude - function to create an ObsPy Magnitude object + :param originID: + :param origintime: + :param cinfo: + :param authority_id: + :return: + ''' + magnresID = createResourceID(origintime, 'mag', authority_id) + magnitude = ope.Magnitude() + magnitude.resource_id = magnresID + magnitude.creation_info = cinfo + magnitude.origin_id = originID + return magnitude + + +def createAmplitude(): + pass + amplID = 'corrampl/' + genID + amplresID = ResourceIdentifier(id=amplID) + amplresID.convertIDToQuakeMLURI(authority_id='BUG') + amplitude = Amplitude() + amplitude.resource_id = amplresID + amplitude.creation_info = self.cinfo + amp = self.data[eventid][phase]['Amplitude'] * 1e-9 + amplitude.generic_amplitude = amp + amplitude.unit = AmplitudeUnit('m/s') + amplitude.magnitude_hint = 'Ml' + amplitude.type = AmplitudeCategory('point') + amplitude.pick_id = pickresID + def getOwner(fn): return pwd.getpwuid(os.stat(fn).st_uid).pw_name