Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/suds.py: 89%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Reader for `SUDS <https://seiscode.iris.washington.edu/projects/suds>`_ 

8waveforms and station information. 

9''' 

10 

11import sys 

12import struct 

13import logging 

14import numpy as num 

15from collections import namedtuple, defaultdict 

16 

17from pyrocko import util, trace, model 

18from .io_common import FileLoadError 

19 

20 

21g_suds_zero = None 

22 

23 

24def get_suds_zero(): 

25 global g_suds_zero 

26 if g_suds_zero is None: 

27 g_suds_zero = util.str_to_time('1970-01-01 00:00:00') 

28 

29 return g_suds_zero 

30 

31 

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

33 

34 

35class SudsError(Exception): 

36 pass 

37 

38 

39class SudsStructBase(object): 

40 

41 @classmethod 

42 def unpack(cls, s): 

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

44 

45 

46class SudsStructtag(SudsStructBase, namedtuple( 

47 'SudsStructtag', 

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

49 

50 __slots__ = () 

51 fmt = '<cchll' 

52 

53 

54class SudsStatident(SudsStructBase, namedtuple( 

55 'SudsStatident', 

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

57 

58 def nslc(self): 

59 return ( 

60 str(self.network.rstrip(b'\0 ').decode('ascii')), 

61 str(self.st_name.rstrip(b'\0 ').decode('ascii')), 

62 '', 

63 str(self.component.rstrip(b'\0 ').decode('ascii'))) 

64 

65 __slots__ = () 

66 fmt = '<4s5sch' 

67 

68 

69class SudsStationcomp(SudsStructBase, namedtuple( 

70 'SudsStationcomp', 

71 'sc_name, azim, incid, st_lat, st_long, elev, enclosure, annotation, ' 

72 'recorder, rockclass, rocktype, sitecondition, sensor_type, ' 

73 'data_type, data_units, polarity, st_status, max_gain, clip_value, ' 

74 'con_mvolts, channel, atod_gain, effective, clock_correct, ' 

75 'station_delay')): 

76 

77 __slots__ = () 

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

79 

80 @classmethod 

81 def unpack(cls, s): 

82 v = struct.unpack(cls.fmt, s) 

83 v = (SudsStatident.unpack(v[0]),) + v[1:] 

84 return cls._make(v) 

85 

86 def to_station(self): 

87 net, sta, loc, cha = self.sc_name.nslc() 

88 station = model.Station( 

89 network=net, 

90 station=sta, 

91 location=loc, 

92 lat=self.st_lat, 

93 lon=self.st_long, 

94 elevation=self.elev) 

95 

96 station.add_channel( 

97 model.Channel( 

98 name=cha, 

99 azimuth=self.azim, 

100 dip=self.incid - 90.)) 

101 

102 return station 

103 

104 

105class SudsDescriptrace(SudsStructBase, namedtuple( 

106 'SudsDescriptrace', 

107 'dt_name, begintime, localtime, datatype, descriptor, digi_by, ' 

108 'processed, length, rate, mindata, maxdata, avenoise, numclip, ' 

109 'time_correct, rate_correct')): 

110 

111 __slots__ = () 

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

113 

114 @classmethod 

115 def unpack(cls, s): 

116 v = struct.unpack(cls.fmt, s) 

117 v = (SudsStatident.unpack(v[0]),) + v[1:] 

118 return cls._make(v) 

119 

120 def to_trace(self, data): 

121 tmin = self.begintime - get_suds_zero() 

122 deltat = 1.0 / self.rate 

123 

124 if data is None: 

125 tmax = tmin + (self.length - 1) * deltat 

126 arr = None 

127 else: 

128 tmax = None 

129 if self.datatype == b'l': 

130 arr = num.frombuffer(data, dtype=num.int32) 

131 elif self.datatype == b'i': 

132 arr = num.frombuffer(data, dtype=num.int16) 

133 elif self.datatype == b'f': 

134 arr = num.frombuffer(data, dtype=num.float32) 

135 elif self.datatype == b'd': 

136 arr = num.frombuffer(data, dtype=num.float64) 

137 else: 

138 raise SudsError( 

139 'data type "%s" not implemented yet' % self.datatype) 

140 

141 if self.length != arr.size: 

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

143 

144 return trace.Trace( 

145 self.dt_name.network.rstrip(b'\0 '), 

146 self.dt_name.st_name.rstrip(b'\0 '), 

147 '', 

148 self.dt_name.component.rstrip(b'\0 '), 

149 ydata=arr, 

150 deltat=deltat, 

151 tmin=tmin, 

152 tmax=tmax) 

153 

154 

155struct_names = { 

156 0: 'no_struct', 

157 1: 'statident', 

158 2: 'structtag', 

159 3: 'terminator', 

160 4: 'equipment', 

161 5: 'stationcomp', 

162 6: 'muxdata', 

163 7: 'descriptrace', 

164 8: 'loctrace', 

165 9: 'calibration', 

166 10: 'feature', 

167 11: 'residual', 

168 12: 'event', 

169 13: 'ev_descript', 

170 14: 'origin', 

171 15: 'error', 

172 16: 'focalmech', 

173 17: 'moment', 

174 18: 'velmodel', 

175 19: 'layers', 

176 20: 'comment', 

177 21: 'profile', 

178 22: 'shotgather', 

179 23: 'calib', 

180 24: 'complex', 

181 25: 'triggers', 

182 26: 'trigsetting', 

183 27: 'eventsetting', 

184 28: 'detector', 

185 29: 'atodinfo', 

186 30: 'timecorrection', 

187 31: 'instrument', 

188 32: 'chanset'} 

189 

190max_struct_type = max(struct_names.keys()) 

191 

192 

193struct_classes = {} 

194for struct_id, struct_name in struct_names.items(): 

195 g = globals() 

196 if struct_id > 0: 

197 class_name = 'Suds' + struct_name.capitalize() 

198 if class_name in g: 

199 struct_classes[struct_id] = g[class_name] 

200 

201 

202def read_suds_struct(f, cls, end_ok=False): 

203 size = struct.calcsize(cls.fmt) 

204 s = f.read(size) 

205 if end_ok and len(s) == 0: 

206 return None 

207 

208 if size != len(s): 

209 raise SudsError('premature end of file') 

210 o = cls.unpack(s) 

211 

212 return o 

213 

214 

215def _iload(filename, load_data=True, want=('traces', 'stations')): 

216 try: 

217 f = open(filename, 'rb') 

218 while True: 

219 tag = read_suds_struct(f, SudsStructtag, end_ok=True) 

220 if tag is None: 

221 break 

222 

223 if tag.struct_type in struct_classes: 

224 cls = struct_classes[tag.struct_type] 

225 if tag.struct_length != struct.calcsize(cls.fmt): 

226 raise SudsError( 

227 'expected and reported struct lengths differ') 

228 

229 s = read_suds_struct(f, cls) 

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

231 station = s.to_station() 

232 yield station 

233 

234 if tag.data_length > 0: 

235 if isinstance(s, SudsDescriptrace) and 'traces' in want: 

236 if load_data: 

237 data = f.read(tag.data_length) 

238 if tag.data_length != len(data): 

239 raise SudsError('premature end of file') 

240 

241 tr = s.to_trace(data) 

242 else: 

243 f.seek(tag.data_length, 1) 

244 tr = s.to_trace(None) 

245 

246 yield tr 

247 else: 

248 f.seek(tag.data_length, 1) 

249 

250 else: 

251 logger.warning( 

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

253 tag.struct_type, 

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

255 

256 f.seek(tag.struct_length, 1) 

257 

258 if tag.data_length > 0: 

259 f.seek(tag.data_length, 1) 

260 

261 except (OSError, SudsError) as e: 

262 raise FileLoadError(e) 

263 

264 finally: 

265 f.close() 

266 

267 

268def iload(filename, load_data=True): 

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

270 yield tr 

271 

272 

273def load_stations(filename): 

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

275 

276 gathered = defaultdict(list) 

277 for s in stations: 

278 nsl = s.nsl() 

279 gathered[nsl].append(s) 

280 

281 stations_out = [] 

282 for nsl, group in gathered.items(): 

283 meta = [ 

284 ('.'.join(s.nsl()) + '.' + s.get_channels()[0].name, 

285 s.lat, s.lon, s.elevation) 

286 for s in group] 

287 

288 util.consistency_check(meta) 

289 

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

291 group[0].set_channels(channels) 

292 stations_out.append(group[0]) 

293 

294 return stations_out 

295 

296 

297def detect(first512): 

298 

299 s = first512[:12] 

300 if len(s) != 12: 

301 return False 

302 

303 tag = SudsStructtag.unpack(s) 

304 if tag.sync != b'S' \ 

305 or tag.machine != b'6' \ 

306 or tag.struct_type < 0 \ 

307 or tag.struct_type > max_struct_type: 

308 

309 return False 

310 

311 return True 

312 

313 

314if __name__ == '__main__': 

315 util.setup_logging('pyrocko.suds') 

316 

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

318 

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

320 

321 for station in stations: 

322 print(station) 

323 

324 trace.snuffle(trs, stations=stations)