|
# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
'''
Module to read and write GSE2.0, GSE2.1, and IMS1.0 files.
'''
from __future__ import print_function, absolute_import
import sys
import re
import logging
from . import util
from .io_common import FileLoadError, FileSaveError
from pyrocko.guts import (
Object, String, StringChoice, Timestamp, Int, Float, List, Bool, Complex,
ValidationError)
try:
range = xrange
except NameError:
pass
logger = logging.getLogger('pyrocko.io.ims')
km = 1000.
nm_per_s = 1.0e-9
g_versions = ('GSE2.0', 'GSE2.1', 'IMS1.0')
g_dialects = ('NOR_NDC', 'USA_DMC')
class SerializeError(Exception):
'''
Raised when serialization of an IMS/GSE2 object fails.
'''
pass
class DeserializeError(Exception):
'''
Raised when deserialization of an IMS/GSE2 object fails.
'''
def __init__(self, *args, **kwargs):
Exception.__init__(self, *args)
self._line_number = None
self._line = None
self._position = kwargs.get('position', None)
self._format = kwargs.get('format', None)
self._version_dialect = None
def set_context(self, line_number, line, version_dialect):
self._line_number = line_number
self._line = line
self._version_dialect = version_dialect
def __str__(self):
lst = [Exception.__str__(self)]
if self._version_dialect is not None:
lst.append('format version: %s' % self._version_dialect[0])
lst.append('dialect: %s' % self._version_dialect[1])
if self._line_number is not None:
lst.append('line number: %i' % self._line_number)
if self._line is not None:
lst.append('line content:\n%s' % (self._line.decode('ascii') or
'*** line is empty ***'))
if self._position is not None:
if self._position[1] is None:
length = max(1, len(self._line or '') - self._position[0])
else:
length = self._position[1]
lst.append(' ' * self._position[0] + '^' * length)
if self._format is not None:
i = 0
f = []
j = 1
for element in self._format:
if element.length != 0:
f.append(' ' * (element.position - i))
if element.length is not None:
f.append(str(j % 10) * element.length)
i = element.position + element.length
else:
f.append(str(j % 10) + '...')
i = element.position + 4
j += 1
lst.append(''.join(f))
return '\n'.join(lst)
def float_or_none(x):
if x.strip():
return float(x)
else:
return None
def int_or_none(x):
if x.strip():
return int(x)
else:
return None
def float_to_string(fmt):
ef = fmt[0]
assert ef in 'ef'
ln, d = map(int, fmt[1:].split('.'))
pfmts = ['%%%i.%i%s' % (ln, dsub, ef) for dsub in range(d, -1, -1)]
blank = b' ' * ln
def func(v):
if v is None:
return blank
for pfmt in pfmts:
s = pfmt % v
if len(s) == ln:
return s.encode('ascii')
raise SerializeError('format="%s", value=%s' % (pfmt, repr(v)))
return func
def int_to_string(fmt):
assert fmt[0] == 'i'
pfmt = '%'+fmt[1:]+'i'
ln = int(fmt[1:])
blank = b' ' * ln
def func(v):
if v is None:
return blank
s = pfmt % v
if len(s) == ln:
return s.encode('ascii')
else:
raise SerializeError('format="%s", value=%s' % (pfmt, repr(v)))
return func
def deserialize_string(fmt):
if fmt.endswith('?'):
def func(s):
if s.strip():
return str(s.rstrip().decode('ascii'))
else:
return None
else:
def func(s):
return str(s.rstrip().decode('ascii'))
return func
def serialize_string(fmt):
if fmt.endswith('+'):
more_ok = True
else:
more_ok = False
fmt = fmt.rstrip('?+')
assert fmt[0] == 'a'
ln = int(fmt[1:])
def func(v):
if v is None:
v = b''
else:
v = v.encode('ascii')
s = v.ljust(ln)
if more_ok or len(s) == ln:
return s
else:
raise SerializeError('max string length: %i, value="%s"' % ln, v)
return func
def rstrip_string(v):
return v.rstrip()
def x_fixed(expect):
def func():
def parse(s):
if s != expect:
raise DeserializeError(
'expected="%s", value="%s"' % (expect, s))
return s
def string(s):
return expect
return parse, string
func.width = len(expect)
func.help_type = 'Keyword: %s' % expect
return func
def x_scaled(fmt, factor):
def func():
to_string = float_to_string(fmt)
def parse(s):
x = float_or_none(s)
if x is None:
return None
else:
return x * factor
def string(v):
if v is None:
return to_string(None)
else:
return to_string(v/factor)
return parse, string
func.width = int(fmt[1:].split('.')[0])
func.help_type = 'float'
return func
def x_int_angle():
def string(v):
if v is None:
return b' '
else:
return ('%3i' % (int(round(v)) % 360)).encode('ascii')
return float_or_none, string
x_int_angle.width = 3
x_int_angle.help_type = 'int [0, 360]'
def x_substitute(value):
def func():
def parse(s):
assert s == b''
return value
def string(s):
return b''
return parse, string
func.width = 0
func.help_type = 'Not present in this file version.'
return func
def fillup_zeros(s, fmt):
s = s.rstrip()
if not s:
return s
if fmt == '%Y/%m/%d %H:%M:%S.3FRAC':
return s + '0000/01/01 00:00:00.000'[len(s):]
elif fmt == '%Y/%m/%d %H:%M:%S.2FRAC':
return s + '0000/01/01 00:00:00.00'[len(s):]
return s
def x_date_time(fmt='%Y/%m/%d %H:%M:%S.3FRAC'):
def parse(s):
s = str(s.decode('ascii'))
try:
s = fillup_zeros(s, fmt)
return util.str_to_time(s, format=fmt)
except Exception:
# iris sets this dummy end dates and they don't fit into 32bit
# time stamps
if fmt[:2] == '%Y' and s[:4] in ('2599', '2045'):
return None
elif fmt[6:8] == '%Y' and s[6:10] in ('2599', '2045'):
return None
raise DeserializeError('expected date, value="%s"' % s)
def string(s):
return util.time_to_str(s, format=fmt).encode('ascii')
return parse, string
x_date_time.width = 23
x_date_time.help_type = 'YYYY/MM/DD HH:MM:SS.FFF'
def x_date():
return x_date_time(fmt='%Y/%m/%d')
x_date.width = 10
x_date.help_type = 'YYYY/MM/DD'
def x_date_iris():
return x_date_time(fmt='%m/%d/%Y')
x_date_iris.width = 10
x_date_iris.help_type = 'MM/DD/YYYY'
def x_date_time_no_seconds():
return x_date_time(fmt='%Y/%m/%d %H:%M')
x_date_time_no_seconds.width = 16
x_date_time_no_seconds.help_type = 'YYYY/MM/DD HH:MM'
def x_date_time_2frac():
return x_date_time(fmt='%Y/%m/%d %H:%M:%S.2FRAC')
x_date_time_2frac.width = 22
x_date_time_2frac.help_type = 'YYYY/MM/DD HH:MM:SS.FF'
def x_yesno():
def parse(s):
if s == b'y':
return True
elif s == b'n':
return False
else:
raise DeserializeError('"y" on "n" expected')
def string(b):
return [b'n', b'y'][int(b)]
return parse, string
x_yesno.width = 1
x_yesno.help_type = 'yes/no'
def optional(x_func):
def func():
parse, string = x_func()
def parse_optional(s):
if s.strip():
return parse(s)
else:
return None
def string_optional(s):
if s is None:
return b' ' * x_func.width
else:
return string(s)
return parse_optional, string_optional
func.width = x_func.width
func.help_type = 'optional %s' % x_func.help_type
return func
class E(object):
def __init__(self, begin, end, fmt, dummy=False):
self.advance = 1
if dummy:
self.advance = 0
self.position = begin - 1
if end is not None:
self.length = end - begin + 1
else:
self.length = None
self.end = end
if isinstance(fmt, str):
t = fmt[0]
if t in 'ef':
self.parse = float_or_none
self.string = float_to_string(fmt)
ln = int(fmt[1:].split('.')[0])
self.help_type = 'float'
elif t == 'a':
self.parse = deserialize_string(fmt)
self.string = serialize_string(fmt)
ln = int(fmt[1:].rstrip('+?'))
self.help_type = 'string'
elif t == 'i':
self.parse = int_or_none
self.string = int_to_string(fmt)
ln = int(fmt[1:])
self.help_type = 'integer'
else:
assert False, 'invalid format: %s' % t
assert self.length is None or ln == self.length, \
'inconsistent length for pos=%i, fmt=%s' \
% (self.position, fmt)
else:
self.parse, self.string = fmt()
self.help_type = fmt.help_type
def end_section(line, extra=None):
if line is None:
return True
ul = line.upper()
return ul.startswith(b'DATA_TYPE') or ul.startswith(b'STOP') or \
(extra is not None and ul.startswith(extra))
class Section(Object):
'''
Base class for top level sections in IMS/GSE2 files.
Sections as understood by this implementation typically correspond to a
DATA_TYPE section in IMS/GSE2 but for some types a finer granularity has
been chosen.
'''
handlers = {} # filled after section have been defined below
@classmethod
def read(cls, reader):
datatype = DataType.read(reader)
reader.pushback()
return Section.handlers[
datatype.type.upper().encode('ascii')].read(reader)
def write_datatype(self, writer):
datatype = DataType(
type=self.keyword.decode('ascii'),
format=writer.version_dialect[0])
datatype.write(writer)
@classmethod
def read_table(cls, reader, expected_header, block_cls, end=end_section):
header = reader.readline()
if not header.upper().startswith(expected_header.upper()):
raise DeserializeError(
'invalid table header line, expected:\n'
'%s\nfound: %s ' % (expected_header, header))
while True:
line = reader.readline()
reader.pushback()
if end(line):
break
yield block_cls.read(reader)
def write_table(self, writer, header, blocks):
writer.writeline(header)
for block in blocks:
block.write(writer)
def get_versioned(x, version_dialect):
if isinstance(x, dict):
for v in (tuple(version_dialect), version_dialect[0], None):
if v in x:
return x[v]
else:
return x
class Block(Object):
'''
Base class for IMS/GSE2 data blocks / lines.
Blocks as understood by this implementation usually correspond to
individual logical lines in the IMS/GSE2 file.
'''
def values(self):
return list(self.T.ivals(self))
@classmethod
def format(cls, version_dialect):
return get_versioned(cls._format, version_dialect)
def serialize(self, version_dialect):
ivalue = 0
out = []
values = self.values()
for element in self.format(version_dialect):
if element.length != 0:
out.append((element.position, element.string(values[ivalue])))
ivalue += element.advance
out.sort()
i = 0
slist = []
for (position, s) in out:
slist.append(b' ' * (position - i))
slist.append(s)
i = position + len(s)
return b''.join(slist)
@classmethod
def deserialize_values(cls, line, version_dialect):
values = []
for element in cls.format(version_dialect):
try:
val = element.parse(
line[element.position:element.end])
if element.advance != 0:
values.append(val)
except Exception:
raise DeserializeError(
'Cannot parse %s' % (
element.help_type),
position=(element.position, element.length),
version=version_dialect[0],
dialect=version_dialect[1],
format=cls.format(version_dialect))
return values
@classmethod
def validated(cls, *args, **kwargs):
obj = cls(*args, **kwargs)
try:
obj.validate()
except ValidationError as e:
raise DeserializeError(str(e))
return obj
@classmethod
def regularized(cls, *args, **kwargs):
obj = cls(*args, **kwargs)
try:
obj.regularize()
except ValidationError as e:
raise DeserializeError(str(e))
return obj
@classmethod
def deserialize(cls, line, version_dialect):
values = cls.deserialize_values(line, version_dialect)
return cls.validated(**dict(zip(cls.T.propnames, values)))
@classmethod
def read(cls, reader):
line = reader.readline()
return cls.deserialize(line, reader.version_dialect)
def write(self, writer):
s = self.serialize(writer.version_dialect)
writer.writeline(s)
class FreeFormatLine(Block):
'''
Base class for IMS/GSE2 free format lines.
'''
@classmethod
def deserialize_values(cls, line, version_dialect):
format = cls.format(version_dialect)
values = line.split(None, len(format)-1)
values_weeded = []
for x, v in zip(format, values):
if isinstance(x, bytes):
if v.upper() != x:
raise DeserializeError(
'expected keyword: %s, found %s' % (x, v.upper()))
else:
if isinstance(x, tuple):
x, (parse, _) = x
v = parse(v)
values_weeded.append((x, v))
values_weeded.sort()
return [str(xv[1].decode('ascii')) for xv in values_weeded]
@classmethod
def deserialize(cls, line, version_dialect):
values = cls.deserialize_values(line, version_dialect)
propnames = cls.T.propnames
stuff = dict(zip(propnames, values))
return cls.regularized(**stuff)
def serialize(self, version_dialect):
names = self.T.propnames
props = self.T.properties
out = []
for x in self.format(version_dialect):
if isinstance(x, bytes):
out.append(x.decode('ascii'))
else:
if isinstance(x, tuple):
x, (_, string) = x
v = string(getattr(self, names[x-1]))
else:
v = getattr(self, names[x-1])
if v is None:
break
out.append(props[x-1].to_save(v))
return ' '.join(out).encode('ascii')
class DataType(Block):
'''
Representation of a DATA_TYPE line.
'''
type = String.T()
subtype = String.T(optional=True)
format = String.T()
subformat = String.T(optional=True)
@classmethod
def deserialize(cls, line, version_dialect):
pat = br'DATA_TYPE +([^ :]+)(:([^ :]+))?( +([^ :]+)(:([^ :]+))?)?'
m = re.match(pat, line)
if not m:
raise DeserializeError('invalid DATA_TYPE line')
return cls.validated(
type=str((m.group(1) or b'').decode('ascii')),
subtype=str((m.group(3) or b'').decode('ascii')),
format=str((m.group(5) or b'').decode('ascii')),
subformat=str((m.group(7) or b'').decode('ascii')))
def serialize(self, version_dialect):
s = self.type
if self.subtype:
s += ':' + self.subtype
f = self.format
if self.subformat:
f += ':' + self.subformat
return ('DATA_TYPE %s %s' % (s, f)).encode('ascii')
@classmethod
def read(cls, reader):
line = reader.readline()
datatype = cls.deserialize(line, reader.version_dialect)
reader.version_dialect[0] = datatype.format
return datatype
def write(self, writer):
s = self.serialize(writer.version_dialect)
writer.version_dialect[0] = self.format
writer.writeline(s)
class FTPFile(FreeFormatLine):
'''
Representation of an FTP_FILE line.
'''
_format = [b'FTP_FILE', 1, 2, 3, 4]
net_address = String.T()
login_mode = StringChoice.T(choices=('USER', 'GUEST'), ignore_case=True)
directory = String.T()
file = String.T()
class WaveformSubformat(StringChoice):
choices = ['INT', 'CM6', 'CM8', 'AU6', 'AU8']
ignore_case = True
class WID2(Block):
'''
Representation of a WID2 line.
'''
_format = [
E(1, 4, x_fixed(b'WID2'), dummy=True),
E(6, 28, x_date_time),
E(30, 34, 'a5'),
E(36, 38, 'a3'),
E(40, 43, 'a4'),
E(45, 47, 'a3'),
E(49, 56, 'i8'),
E(58, 68, 'f11.6'),
E(70, 79, x_scaled('e10.2', nm_per_s)),
E(81, 87, 'f7.3'),
E(89, 94, 'a6?'),
E(96, 100, 'f5.1'),
E(102, 105, 'f4.1')
]
time = Timestamp.T()
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
sub_format = WaveformSubformat.T(default='CM6')
nsamples = Int.T(default=0)
sample_rate = Float.T(default=1.0)
calibration_factor = Float.T(
optional=True,
help='system sensitivity (m/count) at reference period '
'(calibration_period)')
calibration_period = Float.T(
optional=True,
help='calibration reference period [s]')
instrument_type = String.T(
default='', optional=True, help='instrument type (6 characters)')
horizontal_angle = Float.T(
optional=True,
help='horizontal orientation of sensor, clockwise from north [deg]')
vertical_angle = Float.T(
optional=True,
help='vertical orientation of sensor from vertical [deg]')
class OUT2(Block):
'''
Representation of an OUT2 line.
'''
_format = [
E(1, 4, x_fixed(b'OUT2'), dummy=True),
E(6, 28, x_date_time),
E(30, 34, 'a5'),
E(36, 38, 'a3'),
E(40, 43, 'a4'),
E(45, 55, 'f11.3')
]
time = Timestamp.T()
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
duration = Float.T()
class DLY2(Block):
'''
Representation of a DLY2 line.
'''
_format = [
E(1, 4, x_fixed(b'DLY2'), dummy=True),
E(6, 28, x_date_time),
E(30, 34, 'a5'),
E(36, 38, 'a3'),
E(40, 43, 'a4'),
E(45, 55, 'f11.3')
]
time = Timestamp.T()
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
queue_duration = Float.T(help='duration of queue [s]')
class DAT2(Block):
'''
Representation of a DAT2 line.
'''
_format = [
E(1, 4, x_fixed(b'DAT2'), dummy=True)
]
raw_data = List.T(String.T())
@classmethod
def read(cls, reader):
line = reader.readline()
dat2 = cls.deserialize(line, reader.version_dialect)
while True:
line = reader.readline()
if line.upper().startswith(b'CHK2 '):
reader.pushback()
break
else:
if reader._load_data:
dat2.raw_data.append(line.strip())
return dat2
def write(self, writer):
Block.write(self, writer)
for line in self.raw_data:
writer.writeline(line)
class STA2(Block):
'''
Representation of a STA2 line.
'''
_format = [
E(1, 4, x_fixed(b'STA2'), dummy=True),
E(6, 14, 'a9'),
E(16, 24, 'f9.5'),
E(26, 35, 'f10.5'),
E(37, 48, 'a12'),
E(50, 54, x_scaled('f5.3', km)),
E(56, 60, x_scaled('f5.3', km))
]
# the standard requires lat, lon, elevation and depth, we define them as
# optional, however
network = String.T(help='network code (9 characters)')
lat = Float.T(optional=True)
lon = Float.T(optional=True)
coordinate_system = String.T(default='WGS-84')
elevation = Float.T(optional=True, help='elevation [m]')
depth = Float.T(optional=True, help='emplacement depth [m]')
class CHK2(Block):
'''
Representation of a CHK2 line.
'''
_format = [
E(1, 4, x_fixed(b'CHK2'), dummy=True),
E(6, 13, 'i8')
]
checksum = Int.T()
class EID2(Block):
'''
Representation of an EID2 line.
'''
_format = [
E(1, 4, x_fixed(b'EID2'), dummy=True),
E(6, 13, 'a8'),
E(15, 23, 'a9'),
]
event_id = String.T(help='event ID (8 characters)')
bulletin_type = String.T(help='bulletin type (9 characters)')
class BEA2(Block):
'''
Representation of a BEA2 line.
'''
_format = [
E(1, 4, x_fixed(b'BEA2'), dummy=True),
E(6, 17, 'a12'),
E(19, 23, 'f5.1'),
E(25, 29, 'f5.1')]
beam_id = String.T(help='beam ID (12 characters)')
azimuth = Float.T()
slowness = Float.T()
class Network(Block):
'''
Representation of an entry in a NETWORK section.
'''
_format = [
E(1, 9, 'a9'),
E(11, None, 'a64+')]
network = String.T(help='network code (9 characters)')
description = String.T(help='description')
class Station(Block):
'''
Representation of an entry in a STATION section.
'''
_format = {
None: [
E(1, 9, 'a9'),
E(11, 15, 'a5'),
E(17, 20, 'a4'),
E(22, 30, 'f9.5'),
E(32, 41, 'f10.5'),
E(43, 54, 'a12'),
E(56, 60, x_scaled('f5.3', km)),
E(62, 71, x_date),
E(73, 82, optional(x_date))
],
'GSE2.0': [
E(0, -1, x_substitute('')),
E(1, 5, 'a5'),
E(7, 10, 'a4'),
E(12, 20, 'f9.5'),
E(22, 31, 'f10.5'),
E(32, 31, x_substitute('WGS-84')),
E(33, 39, x_scaled('f7.3', km)),
E(41, 50, x_date),
E(52, 61, optional(x_date))]}
_format['IMS1.0', 'USA_DMC'] = list(_format[None])
_format['IMS1.0', 'USA_DMC'][-2:] = [
E(62, 71, x_date_iris),
E(73, 82, optional(x_date_iris))]
network = String.T(help='network code (9 characters)')
station = String.T(help='station code (5 characters)')
type = String.T(
help='station type (4 characters) '
'(1C: single component, 3C: three-component, '
'hfa: high frequency array, lpa: long period array)')
lat = Float.T()
lon = Float.T()
coordinate_system = String.T(default='WGS-84')
elevation = Float.T(help='elevation [m]')
tmin = Timestamp.T()
tmax = Timestamp.T(optional=True)
class Channel(Block):
'''
Representation of an entry in a CHANNEL section.
'''
_format = {
None: [
E(1, 9, 'a9'),
E(11, 15, 'a5'),
E(17, 19, 'a3'),
E(21, 24, 'a4'),
E(26, 34, 'f9.5'),
E(36, 45, 'f10.5'),
E(47, 58, 'a12'),
E(60, 64, x_scaled('f5.3', km)),
E(66, 70, x_scaled('f5.3', km)),
E(72, 77, 'f6.1'),
E(79, 83, 'f5.1'),
E(85, 95, 'f11.6'),
E(97, 102, 'a6'),
E(105, 114, x_date),
E(116, 125, optional(x_date))],
'GSE2.0': [
E(0, -1, x_substitute('')),
E(1, 5, 'a5'),
E(7, 9, 'a3'),
E(11, 14, 'a4'),
E(16, 24, 'f9.5'),
E(26, 35, 'f10.5'),
E(32, 31, x_substitute('WGS-84')),
E(37, 43, x_scaled('f7.3', km)),
E(45, 50, x_scaled('f6.3', km)),
E(52, 57, 'f6.1'),
E(59, 63, 'f5.1'),
E(65, 75, 'f11.6'),
E(77, 83, 'a7'),
E(85, 94, x_date),
E(96, 105, optional(x_date))]}
# norsar plays its own game...
_format['GSE2.0', 'NOR_NDC'] = list(_format['GSE2.0'])
_format['GSE2.0', 'NOR_NDC'][-2:] = [
E(85, 100, x_date_time_no_seconds),
E(102, 117, optional(x_date_time_no_seconds))]
# also iris plays its own game...
_format['IMS1.0', 'USA_DMC'] = list(_format[None])
_format['IMS1.0', 'USA_DMC'][-2:] = [
E(105, 114, x_date_iris),
E(116, 125, optional(x_date_iris))]
network = String.T(help='network code (9 characters)')
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
lat = Float.T(optional=True)
lon = Float.T(optional=True)
coordinate_system = String.T(default='WGS-84')
elevation = Float.T(optional=True, help='elevation [m]')
depth = Float.T(optional=True, help='emplacement depth [m]')
horizontal_angle = Float.T(
optional=True,
help='horizontal orientation of sensor, clockwise from north [deg]')
vertical_angle = Float.T(
optional=True,
help='vertical orientation of sensor from vertical [deg]')
sample_rate = Float.T()
instrument_type = String.T(
default='', optional=True, help='instrument type (6 characters)')
tmin = Timestamp.T()
tmax = Timestamp.T(optional=True)
class BeamGroup(Block):
'''
Representation of an entry in a BEAM group table.
'''
_format = [
E(1, 8, 'a8'),
E(10, 14, 'a5'),
E(16, 18, 'a3'),
E(20, 23, 'a4'),
E(25, 27, 'i3'),
E(29, 37, 'f9.5')]
beam_group = String.T(help='beam group (8 characters)')
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
weight = Int.T(
optional=True,
help='weight used for this component when the beam was formed')
delay = Float.T(
optional=True,
help='beam delay for this component [s] '
'(used for meabs formed by non-plane waves)')
class BeamType(StringChoice):
choices = ['inc', 'coh']
ignore_case = True
class FilterType(StringChoice):
choices = ['BP', 'LP', 'HP', 'BR']
ignore_case = True
class BeamParameters(Block):
'''
Representation of an entry in a BEAM parameters table.
'''
_format = [
E(1, 12, 'a12'),
E(14, 21, 'a8'),
E(23, 25, 'a3'),
E(27, 27, x_yesno),
E(29, 33, 'f5.1'),
E(35, 39, 'f5.3'), # standard says f5.1 -999.0 is vertical beam
E(41, 48, 'a8'),
E(50, 55, 'f6.2'),
E(57, 62, 'f6.2'),
E(64, 65, 'i2'),
E(67, 67, x_yesno),
E(69, 70, 'a2'),
E(72, 81, x_date),
E(83, 92, optional(x_date))]
beam_id = String.T()
beam_group = String.T()
type = BeamType.T()
is_rotated = Bool.T(help='rotation flag')
azimuth = Float.T(
help='azimuth used to steer the beam [deg] (clockwise from North)')
slowness = Float.T(
help='slowness used to steer the beam [s/deg]')
phase = String.T(
help='phase used to set the beam slowness for origin-based beams '
'(8 characters)')
filter_fmin = Float.T(
help='low frequency cut-off for the beam filter [Hz]')
filter_fmax = Float.T(
help='high frequency cut-off for the beam filter [Hz]')
filter_order = Int.T(
help='order of the beam filter')
filter_is_zero_phase = Bool.T(
help='flag to indicate zero-phase filtering')
filter_type = FilterType.T(
help='type of filtering')
tmin = Timestamp.T(
help='start date of beam use')
tmax = Timestamp.T(
optional=True,
help='end date of beam use')
class OutageReportPeriod(Block):
'''
Representation of a the report period of an OUTAGE section.
'''
_format = [
E(1, 18, x_fixed(b'Report period from'), dummy=True),
E(20, 42, x_date_time),
E(44, 45, x_fixed(b'to'), dummy=True),
E(47, 69, x_date_time)]
tmin = Timestamp.T()
tmax = Timestamp.T()
class Outage(Block):
'''
Representation of an entry in the OUTAGE section table.
'''
_format = [
E(1, 9, 'a9'),
E(11, 15, 'a5'),
E(17, 19, 'a3'),
E(21, 24, 'a4'),
E(26, 48, x_date_time),
E(50, 72, x_date_time),
E(74, 83, 'f10.3'),
E(85, None, 'a48+')]
network = String.T(help='network code (9 characters)')
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
tmin = Timestamp.T()
tmax = Timestamp.T()
duration = Float.T()
comment = String.T()
class CAL2(Block):
'''
Representation of a CAL2 line.
'''
_format = {
None: [
E(1, 4, x_fixed(b'CAL2'), dummy=True),
E(6, 10, 'a5'),
E(12, 14, 'a3'),
E(16, 19, 'a4'),
E(21, 26, 'a6'),
E(28, 42, x_scaled('e15.8', nm_per_s)), # standard: e15.2
E(44, 50, 'f7.3'),
E(52, 62, 'f11.5'),
E(64, 79, x_date_time_no_seconds),
E(81, 96, optional(x_date_time_no_seconds))],
'GSE2.0': [
E(1, 4, x_fixed(b'CAL2'), dummy=True),
E(6, 10, 'a5'),
E(12, 14, 'a3'),
E(16, 19, 'a4'),
E(21, 26, 'a6'),
E(28, 37, x_scaled('e10.4', nm_per_s)),
E(39, 45, 'f7.3'),
E(47, 56, 'f10.5'),
E(58, 73, x_date_time_no_seconds),
E(75, 90, optional(x_date_time_no_seconds))]}
station = String.T(help='station code (5 characters)')
channel = String.T(help='channel code (3 characters)')
location = String.T(
default='', optional=True,
help='location code (aux_id, 4 characters)')
instrument_type = String.T(
default='', optional=True, help='instrument type (6 characters)')
calibration_factor = Float.T(
help='system sensitivity (m/count) at reference period '
'(calibration_period)')
calibration_period = Float.T(help='calibration reference period [s]')
sample_rate = Float.T(help='system output sample rate [Hz]')
tmin = Timestamp.T(help='effective start date and time')
tmax = Timestamp.T(optional=True, help='effective end date and time')
comments = List.T(String.T(optional=True))
@classmethod
def read(cls, reader):
lstart = reader.current_line_number()
line = reader.readline()
obj = cls.deserialize(line, reader.version_dialect)
while True:
line = reader.readline()
# make sure all comments are read
if line is None or not line.startswith(b' '):
reader.pushback()
break
obj.append_dataline(line, reader.version_dialect)
obj.comments.extend(reader.get_comments_after(lstart))
return obj
def write(self, writer):
s = self.serialize(writer.version_dialect)
writer.writeline(s)
for c in self.comments:
writer.writeline((' (%s)' % c).encode('ascii'))
class Units(StringChoice):
choices = ['V', 'A', 'C']
ignore_case = True
class Stage(Block):
'''
Base class for IMS/GSE2 response stages.
Available response stages are :py:class:`PAZ2`, :py:class:`FAP2`,
:py:class:`GEN2`, :py:class:`DIG2`, and :py:class:`FIR2`.
'''
stage_number = Int.T(help='stage sequence number')
@classmethod
def read(cls, reader):
lstart = reader.current_line_number()
line = reader.readline()
obj = cls.deserialize(line, reader.version_dialect)
while True:
line = reader.readline()
if line is None or not line.startswith(b' '):
reader.pushback()
break
obj.append_dataline(line, reader.version_dialect)
obj.comments.extend(reader.get_comments_after(lstart))
return obj
def write(self, writer):
line = self.serialize(writer.version_dialect)
writer.writeline(line)
self.write_datalines(writer)
for c in self.comments:
writer.writeline((' (%s)' % c).encode('ascii'))
def write_datalines(self, writer):
pass
class PAZ2Data(Block):
'''
Representation of the complex numbers in PAZ2 sections.
'''
_format = [
E(2, 16, 'e15.8'),
E(18, 32, 'e15.8')]
real = Float.T()
imag = Float.T()
class PAZ2(Stage):
'''
Representation of a PAZ2 line.
'''
_format = {
None: [
E(1, 4, x_fixed(b'PAZ2'), dummy=True),
E(6, 7, 'i2'),
E(9, 9, 'a1'),
E(11, 25, 'e15.8'),
E(27, 30, 'i4'),
E(32, 39, 'f8.3'),
E(41, 43, 'i3'),
E(45, 47, 'i3'),
E(49, None, 'a25+')],
('IMS1.0', 'USA_DMC'): [
E(1, 4, x_fixed(b'PAZ2'), dummy=True),
E(6, 7, 'i2'),
E(9, 9, 'a1'),
E(11, 25, 'e15.8'),
E(27, 30, 'i4'),
E(32, 39, 'f8.3'),
E(40, 42, 'i3'),
E(44, 46, 'i3'),
E(48, None, 'a25+')]}
output_units = Units.T(
help='output units code (V=volts, A=amps, C=counts)')
scale_factor = Float.T(help='scale factor [ouput units/input units]')
decimation = Int.T(optional=True, help='decimation')
correction = Float.T(optional=True, help='group correction applied [s]')
npoles = Int.T(help='number of poles')
nzeros = Int.T(help='number of zeros')
description = String.T(default='', optional=True, help='description')
poles = List.T(Complex.T())
zeros = List.T(Complex.T())
comments = List.T(String.T(optional=True))
def append_dataline(self, line, version_dialect):
d = PAZ2Data.deserialize(line, version_dialect)
v = complex(d.real, d.imag)
i = len(self.poles) + len(self.zeros)
if i < self.npoles:
self.poles.append(v)
elif i < self.npoles + self.nzeros:
self.zeros.append(v)
else:
raise DeserializeError(
'more poles and zeros than expected')
def write_datalines(self, writer):
for pole in self.poles:
PAZ2Data(real=pole.real, imag=pole.imag).write(writer)
for zero in self.zeros:
PAZ2Data(real=zero.real, imag=zero.imag).write(writer)
class FAP2Data(Block):
'''
Representation of the data tuples in FAP2 section.
'''
_format = [
E(2, 11, 'f10.5'),
E(13, 27, 'e15.8'),
E(29, 32, 'i4')]
frequency = Float.T()
amplitude = Float.T()
phase = Float.T()
class FAP2(Stage):
'''
Representation of a FAP2 line.
'''
_format = [
E(1, 4, x_fixed(b'FAP2'), dummy=True),
E(6, 7, 'i2'),
E(9, 9, 'a1'),
E(11, 14, 'i4'),
E(16, 23, 'f8.3'),
E(25, 27, 'i3'),
E(29, 53, 'a25')]
output_units = Units.T(
help='output units code (V=volts, A=amps, C=counts)')
decimation = Int.T(optional=True, help='decimation')
correction = Float.T(help='group correction applied [s]')
ntrip = Int.T(help='number of frequency, amplitude, phase triplets')
description = String.T(default='', optional=True, help='description')
frequencies = List.T(Float.T(), help='frequency [Hz]')
amplitudes = List.T(
Float.T(), help='amplitude [input untits/output units]')
phases = List.T(Float.T(), help='phase delay [degrees]')
comments = List.T(String.T(optional=True))
def append_dataline(self, line, version_dialect):
d = FAP2Data.deserialize(line, version_dialect)
self.frequencies.append(d.frequency)
self.amplitudes.append(d.amplitude)
self.phases.append(d.phase)
def write_datalines(self, writer):
for frequency, amplitude, phase in zip(
self.frequencies, self.amplitudes, self.phases):
FAP2Data(
frequency=frequency,
amplitude=amplitude,
phase=phase).write(writer)
class GEN2Data(Block):
'''
Representation of a data tuple in GEN2 section.
'''
_format = [
E(2, 12, 'f11.5'),
E(14, 19, 'f6.3')]
corner = Float.T(help='corner frequency [Hz]')
slope = Float.T(help='slope above corner [dB/decate]')
class GEN2(Stage):
'''
Representation of a GEN2 line.
'''
_format = [
E(1, 4, x_fixed(b'GEN2'), dummy=True),
E(6, 7, 'i2'),
E(9, 9, 'a1'),
E(11, 25, x_scaled('e15.8', nm_per_s)),
E(27, 33, 'f7.3'),
E(35, 38, 'i4'),
E(40, 47, 'f8.3'),
E(49, 51, 'i3'),
E(53, 77, 'a25')]
output_units = Units.T(
help='output units code (V=volts, A=amps, C=counts)')
calibration_factor = Float.T(
help='system sensitivity (m/count) at reference period '
'(calibration_period)')
calibration_period = Float.T(help='calibration reference period [s]')
decimation = Int.T(optional=True, help='decimation')
correction = Float.T(help='group correction applied [s]')
ncorners = Int.T(help='number of corners')
description = String.T(default='', optional=True, help='description')
corners = List.T(Float.T(), help='corner frequencies [Hz]')
slopes = List.T(Float.T(), help='slopes above corners [dB/decade]')
comments = List.T(String.T(optional=True))
def append_dataline(self, line, version_dialect):
d = GEN2Data.deserialize(line, version_dialect)
self.corners.append(d.corner)
self.slopes.append(d.slope)
def write_datalines(self, writer):
for corner, slope in zip(self.corners, self.slopes):
GEN2Data(corner=corner, slope=slope).write(writer)
class DIG2(Stage):
'''
Representation of a DIG2 line.
'''
_format = [
E(1, 4, x_fixed(b'DIG2'), dummy=True),
E(6, 7, 'i2'),
E(9, 23, 'e15.8'),
E(25, 35, 'f11.5'),
E(37, None, 'a25+')]
sensitivity = Float.T(help='sensitivity [counts/input units]')
sample_rate = Float.T(help='digitizer sample rate [Hz]')
description = String.T(default='', optional=True, help='description')
comments = List.T(String.T(optional=True))
class SymmetryFlag(StringChoice):
choices = ['A', 'B', 'C']
ignore_case = True
class FIR2Data(Block):
'''
Representation of a line of coefficients in a FIR2 section.
'''
_format = [
E(2, 16, 'e15.8'),
E(18, 32, 'e15.8'),
E(34, 48, 'e15.8'),
E(50, 64, 'e15.8'),
E(66, 80, 'e15.8')]
factors = List.T(Float.T())
def values(self):
return self.factors + [None]*(5-len(self.factors))
@classmethod
def deserialize(cls, line, version_dialect):
factors = [v for v in cls.deserialize_values(line, version_dialect)
if v is not None]
return cls.validated(factors=factors)
class FIR2(Stage):
'''
Representation of a FIR2 line.
'''
_format = [
E(1, 4, x_fixed(b'FIR2'), dummy=True),
E(6, 7, 'i2'),
E(9, 18, 'e10.2'),
E(20, 23, 'i4'),
E(25, 32, 'f8.3'),
E(34, 34, 'a1'),
E(36, 39, 'i4'),
E(41, None, 'a25+')]
gain = Float.T(help='filter gain (relative factor, not in dB)')
decimation = Int.T(optional=True, help='decimation')
correction = Float.T(help='group correction applied [s]')
symmetry = SymmetryFlag.T(
help='symmetry flag (A=asymmetric, B=symmetric (odd), '
'C=symmetric (even))')
nfactors = Int.T(help='number of factors')
description = String.T(default='', optional=True, help='description')
comments = List.T(String.T(optional=True))
factors = List.T(Float.T())
def append_dataline(self, line, version_dialect):
d = FIR2Data.deserialize(line, version_dialect)
self.factors.extend(d.factors)
def write_datalines(self, writer):
i = 0
while i < len(self.factors):
FIR2Data(factors=self.factors[i:i+5]).write(writer)
i += 5
class Begin(FreeFormatLine):
'''
Representation of a BEGIN line.
'''
_format = [b'BEGIN', 1]
version = String.T(optional=True)
@classmethod
def read(cls, reader):
line = reader.readline()
obj = cls.deserialize(line, reader.version_dialect)
reader.version_dialect[0] = obj.version
return obj
def write(self, writer):
FreeFormatLine.write(self, writer)
writer.version_dialect[0] = self.version
class MessageType(StringChoice):
choices = ['REQUEST', 'DATA', 'SUBSCRIPTION']
ignore_case = True
class MsgType(FreeFormatLine):
_format = [b'MSG_TYPE', 1]
type = MessageType.T()
class MsgID(FreeFormatLine):
'''
Representation of a MSG_ID line.
'''
_format = [b'MSG_ID', 1, 2]
msg_id_string = String.T()
msg_id_source = String.T(optional=True)
@classmethod
def read(cls, reader):
line = reader.readline()
obj = cls.deserialize(line, reader.version_dialect)
if obj.msg_id_source in g_dialects:
reader.version_dialect[1] = obj.msg_id_source
return obj
def write(self, writer):
FreeFormatLine.write(self, writer)
if self.msg_id_source in g_dialects:
writer.version_dialect[1] = self.msg_id_source
class RefID(FreeFormatLine):
'''
Representation of a REF_ID line.
'''
_format = {
None: [b'REF_ID', 1, 2, 'PART', 3, 'OF', 4],
'GSE2.0': [b'REF_ID', 1]}
msg_id_string = String.T()
msg_id_source = String.T(optional=True)
sequence_number = Int.T(optional=True)
total_number = Int.T(optional=True)
def serialize(self, version_dialect):
out = ['REF_ID', self.msg_id_string]
if self.msg_id_source:
out.append(self.msg_id_source)
i = self.sequence_number
n = self.total_number
if i is not None and n is not None:
out.extend(['PART', str(i), 'OF', str(n)])
return ' '.join(out).encode('ascii')
class LogSection(Section):
'''
Representation of a DATA_TYPE LOG section.
'''
keyword = b'LOG'
lines = List.T(String.T())
@classmethod
def read(cls, reader):
DataType.read(reader)
lines = []
while True:
line = reader.readline()
if end_section(line):
reader.pushback()
break
else:
lines.append(str(line.decode('ascii')))
return cls(lines=lines)
def write(self, writer):
self.write_datatype(writer)
for line in self.lines:
ul = line.upper()
if ul.startswith('DATA_TYPE') or ul.startswith('STOP'):
line = ' ' + line
writer.writeline(str.encode(line))
class ErrorLogSection(LogSection):
'''
Representation of a DATA_TYPE ERROR_LOG section.
'''
keyword = b'ERROR_LOG'
class FTPLogSection(Section):
'''
Representation of a DATA_TYPE FTP_LOG section.
'''
keyword = b'FTP_LOG'
ftp_file = FTPFile.T()
@classmethod
def read(cls, reader):
DataType.read(reader)
ftp_file = FTPFile.read(reader)
return cls(ftp_file=ftp_file)
def write(self, writer):
self.write_datatype(writer)
self.ftp_file.write(writer)
class WID2Section(Section):
'''
Representation of a WID2/STA2/EID2/BEA2/DAT2/CHK2 group.
'''
wid2 = WID2.T()
sta2 = STA2.T(optional=True)
eid2s = List.T(EID2.T())
bea2 = BEA2.T(optional=True)
dat2 = DAT2.T()
chk2 = CHK2.T()
@classmethod
def read(cls, reader):
blocks = dict(eid2s=[])
expect = [(b'WID2 ', WID2, 1)]
if reader.version_dialect[0] == 'GSE2.0':
# should not be there in GSE2.0, but BGR puts it there
expect.append((b'STA2 ', STA2, 0))
else:
expect.append((b'STA2 ', STA2, 1))
expect.extend([
(b'EID2 ', EID2, 0),
(b'BEA2 ', BEA2, 0),
(b'DAT2', DAT2, 1),
(b'CHK2 ', CHK2, 1)])
for k, handler, required in expect:
line = reader.readline()
reader.pushback()
if line is None:
raise DeserializeError('incomplete waveform section')
if line.upper().startswith(k):
block = handler.read(reader)
if k == b'EID2 ':
blocks['eid2s'].append(block)
else:
blocks[str(k.lower().rstrip().decode('ascii'))] = block
else:
if required:
raise DeserializeError('expected %s block' % k)
else:
continue
return cls(**blocks)
def write(self, writer):
self.wid2.write(writer)
if self.sta2:
self.sta2.write(writer)
for eid2 in self.eid2s:
eid2.write(writer)
if self.bea2:
self.bea2.write(writer)
self.dat2.write(writer)
self.chk2.write(writer)
def pyrocko_trace(self, checksum_error='raise'):
from pyrocko import ims_ext, trace
assert checksum_error in ('raise', 'warn', 'ignore')
raw_data = self.dat2.raw_data
nsamples = self.wid2.nsamples
deltat = 1.0 / self.wid2.sample_rate
tmin = self.wid2.time
if self.sta2:
net = self.sta2.network
else:
net = ''
sta = self.wid2.station
loc = self.wid2.location
cha = self.wid2.channel
if raw_data:
ydata = ims_ext.decode_cm6(b''.join(raw_data), nsamples)
if checksum_error != 'ignore':
if ims_ext.checksum(ydata) != self.chk2.checksum:
mess = 'computed checksum value differs from stored value'
if checksum_error == 'raise':
raise DeserializeError(mess)
elif checksum_error == 'warn':
logger.warning(mess)
tmax = None
else:
tmax = tmin + (nsamples - 1) * deltat
ydata = None
return trace.Trace(
net, sta, loc, cha, tmin=tmin, tmax=tmax,
deltat=deltat,
ydata=ydata)
@classmethod
def from_pyrocko_trace(cls, tr,
lat=None, lon=None, elevation=None, depth=None):
from pyrocko import ims_ext
ydata = tr.get_ydata()
raw_data = ims_ext.encode_cm6(ydata)
return cls(
wid2=WID2(
nsamples=tr.data_len(),
sample_rate=1.0 / tr.deltat,
time=tr.tmin,
station=tr.station,
location=tr.location,
channel=tr.channel),
sta2=STA2(
network=tr.network,
lat=lat,
lon=lon,
elevation=elevation,
depth=depth),
dat2=DAT2(
raw_data=[raw_data[i*80:(i+1)*80]
for i in range((len(raw_data)-1)//80 + 1)]),
chk2=CHK2(
checksum=ims_ext.checksum(ydata)))
class OUT2Section(Section):
'''
Representation of a OUT2/STA2 group.
'''
out2 = OUT2.T()
sta2 = STA2.T()
@classmethod
def read(cls, reader):
out2 = OUT2.read(reader)
line = reader.readline()
reader.pushback()
if line.startswith(b'STA2'):
# the spec sais STA2 is mandatory but in practice, it is not
# always there...
sta2 = STA2.read(reader)
else:
sta2 = None
return cls(out2=out2, sta2=sta2)
def write(self, writer):
self.out2.write(writer)
if self.sta2 is not None:
self.sta2.write(writer)
class DLY2Section(Section):
'''
Representation of a DLY2/STA2 group.
'''
dly2 = DLY2.T()
sta2 = STA2.T()
@classmethod
def read(cls, reader):
dly2 = DLY2.read(reader)
sta2 = STA2.read(reader)
return cls(dly2=dly2, sta2=sta2)
def write(self, writer):
self.dly2.write(writer)
self.sta2.write(writer)
class WaveformSection(Section):
'''
Representation of a DATA_TYPE WAVEFORM line.
Any subsequent WID2/OUT2/DLY2 groups are handled as indepenent sections, so
this type just serves as a dummy to read/write the DATA_TYPE WAVEFORM
header.
'''
keyword = b'WAVEFORM'
datatype = DataType.T()
@classmethod
def read(cls, reader):
datatype = DataType.read(reader)
return cls(datatype=datatype)
def write(self, writer):
self.datatype.write(writer)
class TableSection(Section):
'''
Base class for table style sections.
'''
has_data_type_header = True
@classmethod
def read(cls, reader):
if cls.has_data_type_header:
DataType.read(reader)
ts = cls.table_setup
header = get_versioned(ts['header'], reader.version_dialect)
blocks = list(cls.read_table(
reader, header, ts['cls'], end=ts.get('end', end_section)))
return cls(**{ts['attribute']: blocks})
def write(self, writer):
if self.has_data_type_header:
self.write_datatype(writer)
ts = self.table_setup
header = get_versioned(ts['header'], writer.version_dialect)
self.write_table(writer, header, getattr(self, ts['attribute']))
class NetworkSection(TableSection):
'''
Representation of a DATA_TYPE NETWORK section.
'''
keyword = b'NETWORK'
table_setup = dict(
header=b'Net Description',
attribute='networks',
cls=Network)
networks = List.T(Network.T())
class StationSection(TableSection):
'''
Representation of a DATA_TYPE STATION section.
'''
keyword = b'STATION'
table_setup = dict(
header={
None: (
b'Net Sta Type Latitude Longitude Coord '
b'Sys Elev On Date Off Date'),
'GSE2.0': (
b'Sta Type Latitude Longitude Elev On Date '
b'Off Date')},
attribute='stations',
cls=Station)
stations = List.T(Station.T())
class ChannelSection(TableSection):
'''
Representation of a DATA_TYPE CHANNEL section.
'''
keyword = b'CHANNEL'
table_setup = dict(
header={
None: (
b'Net Sta Chan Aux Latitude Longitude Coord Sys'
b' Elev Depth Hang Vang Sample Rate Inst '
b'On Date Off Date'),
'GSE2.0': (
b'Sta Chan Aux Latitude Longitude '
b'Elev Depth Hang Vang Sample_Rate Inst '
b'On Date Off Date')},
attribute='channels',
cls=Channel)
channels = List.T(Channel.T())
class BeamSection(Section):
'''
Representation of a DATA_TYPE BEAM section.
'''
keyword = b'BEAM'
beam_group_header = b'Bgroup Sta Chan Aux Wgt Delay'
beam_parameters_header = b'BeamID Bgroup Btype R Azim Slow '\
b'Phase Flo Fhi O Z F '\
b'On Date Off Date'
group = List.T(BeamGroup.T())
parameters = List.T(BeamParameters.T())
@classmethod
def read(cls, reader):
DataType.read(reader)
def end(line):
return line.upper().startswith(b'BEAMID')
group = list(cls.read_table(reader, cls.beam_group_header, BeamGroup,
end))
parameters = list(cls.read_table(reader, cls.beam_parameters_header,
BeamParameters))
return cls(group=group, parameters=parameters)
def write(self, writer):
self.write_datatype(writer)
self.write_table(writer, self.beam_group_header, self.group)
writer.writeline(b'')
self.write_table(writer, self.beam_parameters_header, self.parameters)
class CAL2Section(Section):
'''
Representation of a CAL2 + stages group in a response section.
'''
cal2 = CAL2.T()
stages = List.T(Stage.T())
@classmethod
def read(cls, reader):
cal2 = CAL2.read(reader)
stages = []
handlers = {
b'PAZ2': PAZ2,
b'FAP2': FAP2,
b'GEN2': GEN2,
b'DIG2': DIG2,
b'FIR2': FIR2}
while True:
line = reader.readline()
reader.pushback()
if end_section(line, b'CAL2'):
break
k = line[:4].upper()
if k in handlers:
stages.append(handlers[k].read(reader))
else:
raise DeserializeError('unexpected line')
return cls(cal2=cal2, stages=stages)
def write(self, writer):
self.cal2.write(writer)
for stage in self.stages:
stage.write(writer)
class ResponseSection(Section):
'''
Representation of a DATA_TYPE RESPONSE line.
Any subsequent CAL2+stages groups are handled as indepenent sections, so
this type just serves as a dummy to read/write the DATA_TYPE RESPONSE
header.
'''
keyword = b'RESPONSE'
datatype = DataType.T()
@classmethod
def read(cls, reader):
datatype = DataType.read(reader)
return cls(datatype=datatype)
def write(self, writer):
self.datatype.write(writer)
class OutageSection(Section):
'''
Representation of a DATA_TYPE OUTAGE section.
'''
keyword = b'OUTAGE'
outages_header = b'NET Sta Chan Aux Start Date Time'\
b' End Date Time Duration Comment'
report_period = OutageReportPeriod.T()
outages = List.T(Outage.T())
@classmethod
def read(cls, reader):
DataType.read(reader)
report_period = OutageReportPeriod.read(reader)
outages = []
outages = list(cls.read_table(reader, cls.outages_header,
Outage))
return cls(
report_period=report_period,
outages=outages)
def write(self, writer):
self.write_datatype(writer)
self.report_period.write(writer)
self.write_table(writer, self.outages_header, self.outages)
class BulletinTitle(Block):
_format = [
E(1, 136, 'a136')]
title = String.T()
g_event_types = dict(
uk='unknown',
ke='known earthquake',
se='suspected earthquake',
kr='known rockburst',
sr='suspected rockburst',
ki='known induced event',
si='suspected induced event',
km='known mine explosion',
sm='suspected mine explosion',
kx='known experimental explosion',
sx='suspected experimental explosion',
kn='known nuclear explosion',
sn='suspected nuclear explosion',
ls='landslide',
de='??',
fe='??',)
class Origin(Block):
_format = [
E(1, 22, x_date_time_2frac),
E(23, 23, 'a1?'),
E(25, 29, 'f5.2'),
E(31, 35, 'f5.2'),
E(37, 44, 'f8.4'),
E(46, 54, 'f9.4'),
E(55, 55, 'a1?'),
E(57, 60, x_scaled('f4.1', km)),
E(62, 66, x_scaled('f5.1', km)),
E(68, 70, x_int_angle),
E(72, 76, x_scaled('f5.1', km)),
E(77, 77, 'a1?'),
E(79, 82, x_scaled('f4.1', km)),
E(84, 87, 'i4'),
E(89, 92, 'i4'),
E(94, 96, x_int_angle),
E(98, 103, 'f6.2'),
E(105, 110, 'f6.2'),
E(112, 112, 'a1?'),
E(114, 114, 'a1?'),
E(116, 117, 'a2?'),
E(119, 127, 'a9'),
E(129, 136, 'a8')]
time = Timestamp.T(
help='epicenter date and time')
time_fixed = StringChoice.T(
choices=['f'],
optional=True,
help='fixed flag, ``"f"`` if fixed origin time solution, '
'``None`` if not')
time_error = Float.T(
optional=True,
help='origin time error [seconds], ``None`` if fixed origin time')
residual = Float.T(
optional=True,
help='root mean square of time residuals [seconds]')
lat = Float.T(
help='latitude')
lon = Float.T(
help='longitude')
lat_lon_fixed = StringChoice.T(
choices=['f'], optional=True,
help='fixed flag, ``"f"`` if fixed epicenter solution, '
'``None`` if not')
ellipse_semi_major_axis = Float.T(
optional=True,
help='semi-major axis of 90% c. i. ellipse or its estimate [m], '
'``None`` if fixed')
ellipse_semi_minor_axis = Float.T(
optional=True,
help='semi-minor axis of 90% c. i. ellipse or its estimate [m], '
'``None`` if fixed')
ellipse_strike = Float.T(
optional=True,
help='strike of 90% c. i. ellipse [0-360], ``None`` if fixed')
depth = Float.T(
help='depth [m]')
depth_fixed = StringChoice.T(
choices=['f', 'd'], optional=True,
help='fixed flag, ``"f"`` fixed depth station, "d" depth phases, '
'``None`` if not fixed depth')
depth_error = Float.T(
optional=True,
help='depth error [m], 90% c. i., ``None`` if fixed depth')
nphases = Int.T(
optional=True,
help='number of defining phases')
nstations = Int.T(
optional=True,
help='number of defining stations')
azimuthal_gap = Float.T(
optional=True,
help='gap in azimuth coverage [deg]')
distance_min = Float.T(
optional=True,
help='distance to closest station [deg]')
distance_max = Float.T(
optional=True,
help='distance to furthest station [deg]')
analysis_type = StringChoice.T(
optional=True,
choices=['a', 'm', 'g'],
help='analysis type, ``"a"`` automatic, ``"m"`` manual, ``"g"`` guess')
location_method = StringChoice.T(
optional=True,
choices=['i', 'p', 'g', 'o'],
help='location method, ``"i"`` inversion, ``"p"`` pattern, '
'``"g"`` ground truth, ``"o"`` other')
event_type = StringChoice.T(
optional=True,
choices=sorted(g_event_types.keys()),
help='event type, ' + ', '.join(
'``"%s"`` %s' % (k, g_event_types[k])
for k in sorted(g_event_types.keys())))
author = String.T(help='author of the origin')
origin_id = String.T(help='origin identification')
class OriginSection(TableSection):
has_data_type_header = False
table_setup = dict(
header={
None: (
b' Date Time Err RMS Latitude Longitude '
b'Smaj Smin Az Depth Err Ndef Nsta Gap mdist Mdist '
b'Qual Author OrigID')},
attribute='origins',
end=lambda line: end_section(line, b'EVENT'),
cls=Origin)
origins = List.T(Origin.T())
class EventTitle(Block):
_format = [
E(1, 5, x_fixed(b'Event'), dummy=True),
E(7, 14, 'a8'),
E(16, 80, 'a65')]
event_id = String.T()
region = String.T()
class EventSection(Section):
'''
Groups Event, Arrival, ...
'''
event_title = EventTitle.T()
origin_section = OriginSection.T()
@classmethod
def read(cls, reader):
event_title = EventTitle.read(reader)
origin_section = OriginSection.read(reader)
return cls(
event_title=event_title,
origin_section=origin_section)
def write(self, writer):
self.event_title.write(writer)
self.origin_section.write(writer)
class EventsSection(Section):
'''
Representation of a DATA_TYPE EVENT section.
'''
keyword = b'EVENT'
bulletin_title = BulletinTitle.T()
event_sections = List.T(EventSection.T())
@classmethod
def read(cls, reader):
DataType.read(reader)
bulletin_title = BulletinTitle.read(reader)
event_sections = []
while True:
line = reader.readline()
reader.pushback()
if end_section(line):
break
if line.upper().startswith(b'EVENT'):
event_sections.append(EventSection.read(reader))
return cls(
bulletin_title=bulletin_title,
event_sections=event_sections,
)
def write(self, writer):
self.write_datatype(writer)
self.bulletin_title.write(writer)
for event_section in self.event_sections:
event_section.write(writer)
class BulletinSection(EventsSection):
'''
Representation of a DATA_TYPE BULLETIN section.
'''
keyword = b'BULLETIN'
for sec in (
LogSection, ErrorLogSection, FTPLogSection, WaveformSection,
NetworkSection, StationSection, ChannelSection, BeamSection,
ResponseSection, OutageSection, EventsSection, BulletinSection):
Section.handlers[sec.keyword] = sec
del sec
class MessageHeader(Section):
'''
Representation of a BEGIN/MSG_TYPE/MSG_ID/REF_ID group.
'''
version = String.T()
type = String.T()
msg_id = MsgID.T(optional=True)
ref_id = RefID.T(optional=True)
@classmethod
def read(cls, reader):
handlers = {
b'BEGIN': Begin,
b'MSG_TYPE': MsgType,
b'MSG_ID': MsgID,
b'REF_ID': RefID}
blocks = {}
while True:
line = reader.readline()
reader.pushback()
ok = False
for k in handlers:
if line.upper().startswith(k):
blocks[k] = handlers[k].read(reader)
ok = True
if not ok:
break
return MessageHeader(
type=blocks[b'MSG_TYPE'].type,
version=blocks[b'BEGIN'].version,
msg_id=blocks.get(b'MSG_ID', None),
ref_id=blocks.get(b'REF_ID', None))
def write(self, writer):
Begin(version=self.version).write(writer)
MsgType(type=self.type).write(writer)
if self.msg_id is not None:
self.msg_id.write(writer)
if self.ref_id is not None:
self.ref_id.write(writer)
def parse_ff_date_time(s):
toks = s.split()
if len(toks) == 2:
sdate, stime = toks
else:
sdate, stime = toks[0], ''
stime += '00:00:00.000'[len(stime):]
return util.str_to_time(
sdate + ' ' + stime, format='%Y/%m/%d %H:%M:%S.3FRAC')
def string_ff_date_time(t):
return util.time_to_str(t, format='%Y/%m/%d %H:%M:%S.3FRAC')
class TimeStamp(FreeFormatLine):
'''
Representation of a TIME_STAMP line.
'''
_format = [b'TIME_STAMP', 1]
value = Timestamp.T()
@classmethod
def deserialize(cls, line, version_dialect):
(s,) = cls.deserialize_values(line, version_dialect)
return cls(value=parse_ff_date_time(s))
def serialize(self, line, version_dialect):
return (
'TIME_STAMP %s' % string_ff_date_time(self.value)).encode('ascii')
class Stop(FreeFormatLine):
'''
Representation of a STOP line.
'''
_format = [b'STOP']
dummy = String.T(optional=True)
class XW01(FreeFormatLine):
'''
Representation of a XW01 line (which is a relict from GSE1).
'''
_format = [b'XW01']
dummy = String.T(optional=True)
re_comment = re.compile(br'^(%(.+)\s*| \((#?)(.+)\)\s*)$')
re_comment_usa_dmc = re.compile(br'^(%(.+)\s*| ?\((#?)(.+)\)\s*)$')
class Reader(object):
def __init__(self, f, load_data=True, version=None, dialect=None):
self._f = f
self._load_data = load_data
self._current_fpos = None
self._current_lpos = None # "physical" line number
self._current_line = None
self._readline_count = 0
self._pushed_back = False
self._handlers = {
b'DATA_TYPE ': Section,
b'WID2 ': WID2Section,
b'OUT2 ': OUT2Section,
b'DLY2 ': DLY2Section,
b'CAL2 ': CAL2Section,
b'BEGIN': MessageHeader,
b'STOP': Stop,
b'XW01': XW01, # for compatibility with BGR dialect
b'HANG:': None, # for compatibility with CNDC
b'VANG:': None,
}
self._comment_lines = []
self._time_stamps = []
self.version_dialect = [version, dialect] # main version, dialect
self._in_garbage = True
def tell(self):
return self._current_fpos
def current_line_number(self):
return self._current_lpos - int(self._pushed_back)
def readline(self):
if self._pushed_back:
self._pushed_back = False
return self._current_line
while True:
self._current_fpos = self._f.tell()
self._current_lpos = self._readline_count + 1
ln = self._f.readline()
self._readline_count += 1
if not ln:
self._current_line = None
return None
lines = [ln.rstrip(b'\n\r')]
while lines[-1].endswith(b'\\'):
lines[-1] = lines[-1][:-1]
ln = self._f.readline()
self._readline_count += 1
lines.append(ln.rstrip(b'\n\r'))
self._current_line = b''.join(lines)
if self.version_dialect[1] == 'USA_DMC':
m_comment = re_comment_usa_dmc.match(self._current_line)
else:
m_comment = re_comment.match(self._current_line)
if not self._current_line.strip():
pass
elif m_comment:
comment_type = None
if m_comment.group(3) == b'#':
comment_type = 'ISF'
elif m_comment.group(4) is not None:
comment_type = 'IMS'
comment = m_comment.group(2) or m_comment.group(4)
self._comment_lines.append(
(self._current_lpos, comment_type,
str(comment.decode('ascii'))))
elif self._current_line[:10].upper() == b'TIME_STAMP':
self._time_stamps.append(
TimeStamp.deserialize(
self._current_line, self.version_dialect))
else:
return self._current_line
def get_comments_after(self, lpos):
comments = []
i = len(self._comment_lines) - 1
while i >= 0:
if self._comment_lines[i][0] <= lpos:
break
comments.append(self._comment_lines[i][-1])
i -= 1
return comments
def pushback(self):
assert not self._pushed_back
self._pushed_back = True
def __iter__(self):
return self
def next(self):
return self.__next__()
def __next__(self):
try:
while True:
line = self.readline()
if line is None:
raise StopIteration()
ignore = False
for k in self._handlers:
if line.upper().startswith(k):
if self._handlers[k] is None:
ignore = True
break
self.pushback()
sec = self._handlers[k].read(self)
if isinstance(sec, Stop):
self._in_garbage = True
else:
self._in_garbage = False
return sec
if not self._in_garbage and not ignore:
raise DeserializeError('unexpected line')
except DeserializeError as e:
e.set_context(
self._current_lpos,
self._current_line,
self.version_dialect)
raise
class Writer(object):
def __init__(self, f, version='GSE2.1', dialect=None):
self._f = f
self.version_dialect = [version, dialect]
def write(self, section):
section.write(self)
def writeline(self, line):
self._f.write(line.rstrip())
self._f.write(b'\n')
def write_string(sections):
from io import BytesIO
f = BytesIO()
w = Writer(f)
for section in sections:
w.write(section)
return f.getvalue()
def iload_fh(f, **kwargs):
'''
Load IMS/GSE2 records from open file handle.
'''
try:
r = Reader(f, **kwargs)
for section in r:
yield section
except DeserializeError as e:
raise FileLoadError(e)
def iload_string(s, **kwargs):
'''
Read IMS/GSE2 sections from bytes string.
'''
from io import BytesIO
f = BytesIO(s)
return iload_fh(f, **kwargs)
iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family(
iload_fh, 'IMS/GSE2', ':py:class:`Section`')
def dump_fh(sections, f):
'''
Dump IMS/GSE2 sections to open file handle.
'''
try:
w = Writer(f)
for section in sections:
w.write(section)
except SerializeError as e:
raise FileSaveError(e)
def dump_string(sections):
'''
Write IMS/GSE2 sections to string.
'''
from io import BytesIO
f = BytesIO()
dump_fh(sections, f)
return f.getvalue()
if __name__ == '__main__':
from optparse import OptionParser
usage = 'python -m pyrocko.ims <filenames>'
util.setup_logging('pyrocko.ims.__main__', 'warning')
description = '''
Read and print IMS/GSE2 records.
'''
parser = OptionParser(
usage=usage,
description=description,
formatter=util.BetterHelpFormatter())
parser.add_option(
'--version',
dest='version',
choices=g_versions,
help='inial guess for version')
parser.add_option(
'--dialect',
dest='dialect',
choices=g_dialects,
help='inial guess for dialect')
parser.add_option(
'--load-data',
dest='load_data',
action='store_true',
help='unpack data samples')
parser.add_option(
'--out-version',
dest='out_version',
choices=g_versions,
help='output format version')
parser.add_option(
'--out-dialect',
dest='out_dialect',
choices=g_dialects,
help='output format dialect')
(options, args) = parser.parse_args(sys.argv[1:])
for fn in args:
with open(fn, 'rb') as f:
r = Reader(f, load_data=options.load_data,
version=options.version, dialect=options.dialect)
w = None
if options.out_version is not None:
w = Writer(
sys.stdout, version=options.out_version,
dialect=options.out_dialect)
for sec in r:
if not w:
print(sec)
else:
w.write(sec)
if isinstance(sec, WID2Section) and options.load_data:
tr = sec.pyrocko_trace(checksum_error='warn')
|