1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import division, absolute_import, print_function
7import sys
8import calendar
9import numpy as num
11from pyrocko import util, trace
12from pyrocko.util import unpack_fixed
13from .io_common import FileLoadError
16class SeisanFileError(Exception):
17 pass
20def read_file_header(f, npad=4):
21 header_infos = []
23 nlines = 12
24 iline = 0
25 while iline < nlines:
26 f.read(npad)
27 d = f.read(80)
28 d = str(d.decode('ascii'))
29 f.read(npad)
31 if iline == 0:
32 net_name, nchannels, ear, doy, mon, day, hr, min, secs, tlen = \
33 unpack_fixed(
34 'x1,a29,i3,i3,x1,i3,x1,i2,x1,i2,x1,i2,x1,i2,x1,f6,x1,f9',
35 d)
37 year = 1900 + ear
38 tmin = util.to_time_float(
39 calendar.timegm((year, mon, day, hr, min, secs)))
40 header_infos.append(
41 (net_name, nchannels, util.time_to_str(tmin)))
43 if nchannels > 30:
44 nlines += (nchannels - 31)//3 + 1
46 if iline >= 2:
47 for j in range(3):
48 s = d[j*26:(j+1)*26]
49 if s.strip():
50 sta1, cha1, cha2, sta2, toffset, tlen = unpack_fixed(
51 'x1,a4,a2,x1,a1,a1,f7,x1,f8', s)
53 sta = sta1 + sta2
54 cha = cha1 + cha2
55 header_infos.append(
56 (sta, cha, toffset, tlen))
58 iline += 1
59 return header_infos
62class EOF(Exception):
63 pass
66def read_channel_header(f, npad=4):
67 x = f.read(npad)
68 if len(x) == 0:
69 raise EOF()
71 d = f.read(1040)
72 d = str(d.decode('ascii'))
73 f.read(npad)
75 sta, cha1, loc1, cha2, ear, loc2, doy, net1, mon, net2, day, hr, min, \
76 tflag, secs, rate, nsamples, lat, lon, elevation, gain_flag, \
77 sample_bytes, response_flag1, response_flag2 = unpack_fixed(
78 'a5,a2,a1,a1,i3,a1,i3,a1,i2,a1,i2,x1,i2,x1,i2,a1,f6,x1,f7,i7,x1,'
79 'f8?,x1,f9?,x1,f5?,a1,i1,a1,a1', d[:79])
81 gain = 1
82 if gain_flag:
83 gain = unpack_fixed('f12', d[147:159])
85 cha = cha1+cha2
86 loc = loc1+loc2
87 net = net1+net2
88 tmin = util.to_time_float(
89 calendar.timegm((1900+ear, mon, day, hr, min, secs)))
90 deltat = 1./rate
92 return (net, sta, loc, cha,
93 tmin, tflag, deltat, nsamples, sample_bytes,
94 lat, lon, elevation, gain)
97def read_channel_data(
98 f, endianness, sample_bytes, nsamples, gain, load_data=True, npad=4):
100 if not load_data:
101 f.seek(sample_bytes*nsamples + 2*npad, 1)
102 return None
104 else:
105 f.read(npad)
106 data = num.fromfile(
107 f,
108 dtype=num.dtype('%si%i' % (endianness, sample_bytes)),
109 count=nsamples).astype('i%i' % sample_bytes)
111 f.read(npad)
112 data *= gain
113 return data
116def iload(filename, load_data=True, subformat='l4'):
118 try:
119 npad = 4
120 if subformat is not None:
121 try:
122 endianness = {'l': '<', 'b': '>'}[subformat[0]]
123 if len(subformat) > 1:
124 npad = int(subformat[1:])
125 except Exception:
126 raise SeisanFileError(
127 'Bad subformat specification: "%s"' % subformat)
128 else:
129 endianness = '<'
131 with open(filename, 'rb') as f:
132 try:
133 read_file_header(f, npad=npad)
135 except util.UnpackError as e:
136 raise SeisanFileError(
137 'Error loading header from file %s: %s'
138 % (filename, str(e)))
140 try:
141 itrace = 0
142 while True:
143 try:
144 (net, sta, loc, cha, tmin, tflag, deltat, nsamples,
145 sample_bytes, lat, lon, elevation, gain) \
146 = read_channel_header(f, npad=npad)
148 data = read_channel_data(
149 f, endianness, sample_bytes, nsamples,
150 gain, load_data,
151 npad=npad)
153 tmax = None
154 if data is None:
155 tmax = tmin + (nsamples-1)*deltat
157 t = trace.Trace(
158 net, sta, loc, cha, tmin=tmin, tmax=tmax,
159 deltat=deltat, ydata=data)
161 yield t
163 except util.UnpackError as e:
164 raise SeisanFileError(
165 'Error loading trace %i from file %s: %s' % (
166 itrace, filename, str(e)))
168 itrace += 1
170 except EOF:
171 pass
173 except (OSError, SeisanFileError) as e:
174 raise FileLoadError(e)
177if __name__ == '__main__':
178 fn = sys.argv[1]
179 for tr in iload(fn):
180 print(tr)