1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5from __future__ import absolute_import, division, print_function 

6 

7import sys 

8import struct 

9import logging 

10import numpy as num 

11from collections import namedtuple, defaultdict 

12 

13from pyrocko import util, trace, model 

14from .io_common import FileLoadError 

15 

16 

17g_suds_zero = None 

18 

19 

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') 

24 

25 return g_suds_zero 

26 

27 

28logger = logging.getLogger('pyrocko.io.suds') 

29 

30 

31class SudsError(Exception): 

32 pass 

33 

34 

35class SudsStructBase(object): 

36 

37 @classmethod 

38 def unpack(cls, s): 

39 return cls._make(struct.unpack(cls.fmt, s)) 

40 

41 

42class SudsStructtag(SudsStructBase, namedtuple( 

43 'SudsStructtag', 

44 'sync, machine, struct_type, struct_length, data_length')): 

45 

46 __slots__ = () 

47 fmt = '<cchll' 

48 

49 

50class SudsStatident(SudsStructBase, namedtuple( 

51 'SudsStatident', 

52 'network, st_name, component, inst_type')): 

53 

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'))) 

60 

61 __slots__ = () 

62 fmt = '<4s5sch' 

63 

64 

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')): 

72 

73 __slots__ = () 

74 fmt = '<%ishhddfcccchccccccfffhhlff' % struct.calcsize(SudsStatident.fmt) 

75 

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) 

81 

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) 

91 

92 station.add_channel( 

93 model.Channel( 

94 name=cha, 

95 azimuth=self.azim, 

96 dip=self.incid - 90.)) 

97 

98 return station 

99 

100 

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')): 

106 

107 __slots__ = () 

108 fmt = '<%isdhcchhlffffldf' % struct.calcsize(SudsStatident.fmt) 

109 

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) 

115 

116 def to_trace(self, data): 

117 tmin = self.begintime - get_suds_zero() 

118 deltat = 1.0 / self.rate 

119 

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) 

136 

137 if self.length != arr.size: 

138 raise SudsError('found and reported number of samples differ') 

139 

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) 

149 

150 

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'} 

185 

186max_struct_type = max(struct_names.keys()) 

187 

188 

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] 

196 

197 

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 

203 

204 if size != len(s): 

205 raise SudsError('premature end of file') 

206 o = cls.unpack(s) 

207 

208 return o 

209 

210 

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 

218 

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') 

224 

225 s = read_suds_struct(f, cls) 

226 if isinstance(s, SudsStationcomp) and 'stations' in want: 

227 station = s.to_station() 

228 yield station 

229 

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') 

236 

237 tr = s.to_trace(data) 

238 else: 

239 f.seek(tag.data_length, 1) 

240 tr = s.to_trace(None) 

241 

242 yield tr 

243 else: 

244 f.seek(tag.data_length, 1) 

245 

246 else: 

247 logger.warning( 

248 'skipping unsupported SUDS struct type %s (%s)' % ( 

249 tag.struct_type, 

250 struct_names.get(tag.struct_type, '?'))) 

251 

252 f.seek(tag.struct_length, 1) 

253 

254 if tag.data_length > 0: 

255 f.seek(tag.data_length, 1) 

256 

257 except (OSError, SudsError) as e: 

258 raise FileLoadError(e) 

259 

260 finally: 

261 f.close() 

262 

263 

264def iload(filename, load_data=True): 

265 for tr in _iload(filename, load_data=load_data, want=('traces',)): 

266 yield tr 

267 

268 

269def load_stations(filename): 

270 stations = list(_iload(filename, load_data=False, want=('stations',))) 

271 

272 gathered = defaultdict(list) 

273 for s in stations: 

274 nsl = s.nsl() 

275 gathered[nsl].append(s) 

276 

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] 

283 

284 util.consistency_check(meta) 

285 

286 channels = [s.get_channels()[0] for s in group] 

287 group[0].set_channels(channels) 

288 stations_out.append(group[0]) 

289 

290 return stations_out 

291 

292 

293def detect(first512): 

294 

295 s = first512[:12] 

296 if len(s) != 12: 

297 return False 

298 

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: 

304 

305 return False 

306 

307 return True 

308 

309 

310if __name__ == '__main__': 

311 util.setup_logging('pyrocko.suds') 

312 

313 trs = list(iload(sys.argv[1], 'rb')) 

314 

315 stations = load_stations(sys.argv[1]) 

316 

317 for station in stations: 

318 print(station) 

319 

320 trace.snuffle(trs, stations=stations)