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