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