1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import sys
8import struct
9import logging
10import numpy as num
11from collections import namedtuple, defaultdict
13from pyrocko import util, trace, model
14from .io_common import FileLoadError
17g_suds_zero = None
20def get_suds_zero():
21 global g_suds_zero
22 if g_suds_zero is None:
23 g_suds_zero = util.str_to_time('1970-01-01 00:00:00')
25 return g_suds_zero
28logger = logging.getLogger('pyrocko.io.suds')
31class SudsError(Exception):
32 pass
35class SudsStructBase(object):
37 @classmethod
38 def unpack(cls, s):
39 return cls._make(struct.unpack(cls.fmt, s))
42class SudsStructtag(SudsStructBase, namedtuple(
43 'SudsStructtag',
44 'sync, machine, struct_type, struct_length, data_length')):
46 __slots__ = ()
47 fmt = '<cchll'
50class SudsStatident(SudsStructBase, namedtuple(
51 'SudsStatident',
52 'network, st_name, component, inst_type')):
54 def nslc(self):
55 return (
56 str(self.network.rstrip(b'\0 ').decode('ascii')),
57 str(self.st_name.rstrip(b'\0 ').decode('ascii')),
58 '',
59 str(self.component.rstrip(b'\0 ').decode('ascii')))
61 __slots__ = ()
62 fmt = '<4s5sch'
65class SudsStationcomp(SudsStructBase, namedtuple(
66 'SudsStationcomp',
67 'sc_name, azim, incid, st_lat, st_long, elev, enclosure, annotation, '
68 'recorder, rockclass, rocktype, sitecondition, sensor_type, '
69 'data_type, data_units, polarity, st_status, max_gain, clip_value, '
70 'con_mvolts, channel, atod_gain, effective, clock_correct, '
71 'station_delay')):
73 __slots__ = ()
74 fmt = '<%ishhddfcccchccccccfffhhlff' % struct.calcsize(SudsStatident.fmt)
76 @classmethod
77 def unpack(cls, s):
78 v = struct.unpack(cls.fmt, s)
79 v = (SudsStatident.unpack(v[0]),) + v[1:]
80 return cls._make(v)
82 def to_station(self):
83 net, sta, loc, cha = self.sc_name.nslc()
84 station = model.Station(
85 network=net,
86 station=sta,
87 location=loc,
88 lat=self.st_lat,
89 lon=self.st_long,
90 elevation=self.elev)
92 station.add_channel(
93 model.Channel(
94 name=cha,
95 azimuth=self.azim,
96 dip=self.incid - 90.))
98 return station
101class SudsDescriptrace(SudsStructBase, namedtuple(
102 'SudsDescriptrace',
103 'dt_name, begintime, localtime, datatype, descriptor, digi_by, '
104 'processed, length, rate, mindata, maxdata, avenoise, numclip, '
105 'time_correct, rate_correct')):
107 __slots__ = ()
108 fmt = '<%isdhcchhlffffldf' % struct.calcsize(SudsStatident.fmt)
110 @classmethod
111 def unpack(cls, s):
112 v = struct.unpack(cls.fmt, s)
113 v = (SudsStatident.unpack(v[0]),) + v[1:]
114 return cls._make(v)
116 def to_trace(self, data):
117 tmin = self.begintime - get_suds_zero()
118 deltat = 1.0 / self.rate
120 if data is None:
121 tmax = tmin + (self.length - 1) * deltat
122 arr = None
123 else:
124 tmax = None
125 if self.datatype == b'l':
126 arr = num.fromstring(data, dtype=num.int32)
127 elif self.datatype == b'i':
128 arr = num.fromstring(data, dtype=num.int16)
129 elif self.datatype == b'f':
130 arr = num.fromstring(data, dtype=num.float32)
131 elif self.datatype == b'd':
132 arr = num.fromstring(data, dtype=num.float64)
133 else:
134 raise SudsError(
135 'data type "%s" not implemented yet' % self.datatype)
137 if self.length != arr.size:
138 raise SudsError('found and reported number of samples differ')
140 return trace.Trace(
141 self.dt_name.network.rstrip(b'\0 '),
142 self.dt_name.st_name.rstrip(b'\0 '),
143 '',
144 self.dt_name.component.rstrip(b'\0 '),
145 ydata=arr,
146 deltat=deltat,
147 tmin=tmin,
148 tmax=tmax)
151struct_names = {
152 0: 'no_struct',
153 1: 'statident',
154 2: 'structtag',
155 3: 'terminator',
156 4: 'equipment',
157 5: 'stationcomp',
158 6: 'muxdata',
159 7: 'descriptrace',
160 8: 'loctrace',
161 9: 'calibration',
162 10: 'feature',
163 11: 'residual',
164 12: 'event',
165 13: 'ev_descript',
166 14: 'origin',
167 15: 'error',
168 16: 'focalmech',
169 17: 'moment',
170 18: 'velmodel',
171 19: 'layers',
172 20: 'comment',
173 21: 'profile',
174 22: 'shotgather',
175 23: 'calib',
176 24: 'complex',
177 25: 'triggers',
178 26: 'trigsetting',
179 27: 'eventsetting',
180 28: 'detector',
181 29: 'atodinfo',
182 30: 'timecorrection',
183 31: 'instrument',
184 32: 'chanset'}
186max_struct_type = max(struct_names.keys())
189struct_classes = {}
190for struct_id, struct_name in struct_names.items():
191 g = globals()
192 if struct_id > 0:
193 class_name = 'Suds' + struct_name.capitalize()
194 if class_name in g:
195 struct_classes[struct_id] = g[class_name]
198def read_suds_struct(f, cls, end_ok=False):
199 size = struct.calcsize(cls.fmt)
200 s = f.read(size)
201 if end_ok and len(s) == 0:
202 return None
204 if size != len(s):
205 raise SudsError('premature end of file')
206 o = cls.unpack(s)
208 return o
211def _iload(filename, load_data=True, want=('traces', 'stations')):
212 try:
213 f = open(filename, 'rb')
214 while True:
215 tag = read_suds_struct(f, SudsStructtag, end_ok=True)
216 if tag is None:
217 break
219 if tag.struct_type in struct_classes:
220 cls = struct_classes[tag.struct_type]
221 if tag.struct_length != struct.calcsize(cls.fmt):
222 raise SudsError(
223 'expected and reported struct lengths differ')
225 s = read_suds_struct(f, cls)
226 if isinstance(s, SudsStationcomp) and 'stations' in want:
227 station = s.to_station()
228 yield station
230 if tag.data_length > 0:
231 if isinstance(s, SudsDescriptrace) and 'traces' in want:
232 if load_data:
233 data = f.read(tag.data_length)
234 if tag.data_length != len(data):
235 raise SudsError('premature end of file')
237 tr = s.to_trace(data)
238 else:
239 f.seek(tag.data_length, 1)
240 tr = s.to_trace(None)
242 yield tr
243 else:
244 f.seek(tag.data_length, 1)
246 else:
247 logger.warning(
248 'skipping unsupported SUDS struct type %s (%s)' % (
249 tag.struct_type,
250 struct_names.get(tag.struct_type, '?')))
252 f.seek(tag.struct_length, 1)
254 if tag.data_length > 0:
255 f.seek(tag.data_length, 1)
257 except (OSError, SudsError) as e:
258 raise FileLoadError(e)
260 finally:
261 f.close()
264def iload(filename, load_data=True):
265 for tr in _iload(filename, load_data=load_data, want=('traces',)):
266 yield tr
269def load_stations(filename):
270 stations = list(_iload(filename, load_data=False, want=('stations',)))
272 gathered = defaultdict(list)
273 for s in stations:
274 nsl = s.nsl()
275 gathered[nsl].append(s)
277 stations_out = []
278 for nsl, group in gathered.items():
279 meta = [
280 ('.'.join(s.nsl()) + '.' + s.get_channels()[0].name,
281 s.lat, s.lon, s.elevation)
282 for s in group]
284 util.consistency_check(meta)
286 channels = [s.get_channels()[0] for s in group]
287 group[0].set_channels(channels)
288 stations_out.append(group[0])
290 return stations_out
293def detect(first512):
295 s = first512[:12]
296 if len(s) != 12:
297 return False
299 tag = SudsStructtag.unpack(s)
300 if tag.sync != b'S' \
301 or tag.machine != b'6' \
302 or tag.struct_type < 0 \
303 or tag.struct_type > max_struct_type:
305 return False
307 return True
310if __name__ == '__main__':
311 util.setup_logging('pyrocko.suds')
313 trs = list(iload(sys.argv[1], 'rb'))
315 stations = load_stations(sys.argv[1])
317 for station in stations:
318 print(station)
320 trace.snuffle(trs, stations=stations)