From f6bf37c92035085adbcf6573bb3b408f20508353 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Sat, 7 Feb 2015 09:12:58 +0100 Subject: [PATCH] new package io.py which should contain all import and export routines written by ourselves, such as reading old PILOT phase and location information file; implementation of the readPILOTevent function; new routines in utils: createArrival will be split into two functions: createPick and createArrival; also planned: createOrigin, createAmplitude and createMagnitude as well as giving createEvent functionality --- pylot/RELEASE-VERSION | 2 +- pylot/core/read/io.py | 170 ++++++++++++++++++++++++++++++++++++ pylot/core/util/__init__.py | 2 + pylot/core/util/errors.py | 2 +- pylot/core/util/utils.py | 34 +++++++- 5 files changed, 206 insertions(+), 4 deletions(-) create mode 100644 pylot/core/read/io.py diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION index 52a4577d..05b25a55 100644 --- a/pylot/RELEASE-VERSION +++ b/pylot/RELEASE-VERSION @@ -1 +1 @@ -ef50-dirty +3667-dirty diff --git a/pylot/core/read/io.py b/pylot/core/read/io.py new file mode 100644 index 00000000..ff6ee866 --- /dev/null +++ b/pylot/core/read/io.py @@ -0,0 +1,170 @@ +#!/usr/bin/env python +# -*- 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 + +def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): + """ + readPILOTEvent - function + + Reads Matlab PHASES and LOC files written by Matlab versions of PILOT and + converts the data into an ObsPy Event object which is returned to the + calling program. + + :rtype : ~obspy.core.event.Event + :param eventID: + :param authority: + :param kwargs: + :param phasfn: filename of the old PILOT Matlab PHASES file + :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)) + phasauthor = getOwner(phasfn) + else: + phases = None + phasctime = None + phasauthor = None + if locfn is not None and os.path.isfile(locfn): + loc = sio.loadmat(locfn) + locctime = UTCDateTime(os.path.getmtime(locfn)) + locauthor = getOwner(locfn) + else: + loc = None + locctime = None + locauthor = None + pickcinfo = ope.CreationInfo(agency_id=authority_id, + author=phasauthor, + creation_time=phasctime) + 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('.') + + # retrieve date information + julday = int(idsplit[1]) + year = int(idsplit[2]) + hour = int(loc['hh']) + minute = int(loc['mm']) + second = int(loc['ss']) + + if UTCDateTime(year=year+2000) < UTCDateTime.utcnow(): + year += 2000 + else: + year += 1900 + + eventDate = UTCDateTime(year=year, julday=julday, hour=hour, + minute=minute, second=second) + + 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 + 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])} + 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) + + for picktime, phase in [(ppicktime, 'P'), (spicktime, 'S')]: + if phase == 'P': + 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) + event.picks.append(pick) + + 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) + + 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 + + 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)) + + + diff --git a/pylot/core/util/__init__.py b/pylot/core/util/__init__.py index a6b5d183..4a856b7a 100755 --- a/pylot/core/util/__init__.py +++ b/pylot/core/util/__init__.py @@ -5,6 +5,8 @@ from pylot.core.util.errors import 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 diff --git a/pylot/core/util/errors.py b/pylot/core/util/errors.py index 0b2a4780..e2157135 100644 --- a/pylot/core/util/errors.py +++ b/pylot/core/util/errors.py @@ -10,4 +10,4 @@ class OptionsError(Exception): pass class FormatError(Exception): - pass + pass \ No newline at end of file diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 9748d7eb..1ec04f60 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -2,8 +2,10 @@ # # -*- coding: utf-8 -*- +import os +import pwd import re -from obspy.core.event import * +import obspy.core.event as ope def fnConstructor(s): @@ -19,6 +21,34 @@ def fnConstructor(s): return fn def createEvent(origintime, latitude, longitude, depth, **kwargs): - evt = Event() + 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) + 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:/') + + arriID = 'arrival/' + eventnum + '/' + station + '/{0}'.format(phase) + arriresID = ope.ResourceIdentifier(id=arriID) + arriresID.convertIDToQuakeMLURI(authority_id=authority_id) + arrival = ope.Arrival() + 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'] + +def getOwner(fn): + return pwd.getpwuid(os.stat(fn).st_uid).pw_name +