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
This commit is contained in:
170
pylot/core/read/io.py
Normal file
170
pylot/core/read/io.py
Normal file
@@ -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))
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user