1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6from struct import unpack 

7import os 

8import re 

9import math 

10import logging 

11 

12from pyrocko import trace 

13from pyrocko.util import reuse, ensuredirs 

14from .io_common import FileLoadError, FileSaveError 

15 

16logger = logging.getLogger('pyrocko.io.mseed') 

17 

18MSEED_HEADER_BYTES = 64 

19VALID_RECORD_LENGTHS = tuple(2**exp for exp in range(8, 20)) 

20 

21 

22class CodeTooLong(FileSaveError): 

23 pass 

24 

25 

26def iload(filename, load_data=True, offset=0, segment_size=0, nsegments=0): 

27 from pyrocko import mseed_ext 

28 

29 have_zero_rate_traces = False 

30 try: 

31 isegment = 0 

32 while isegment < nsegments or nsegments == 0: 

33 tr_tuples = mseed_ext.get_traces( 

34 filename, load_data, offset, segment_size) 

35 

36 if not tr_tuples: 

37 break 

38 

39 for tr_tuple in tr_tuples: 

40 network, station, location, channel = tr_tuple[1:5] 

41 tmin = float(tr_tuple[5])/float(mseed_ext.HPTMODULUS) 

42 tmax = float(tr_tuple[6])/float(mseed_ext.HPTMODULUS) 

43 try: 

44 deltat = reuse(1.0/float(tr_tuple[7])) 

45 except ZeroDivisionError: 

46 have_zero_rate_traces = True 

47 continue 

48 

49 ydata = tr_tuple[8] 

50 

51 tr = trace.Trace( 

52 network.strip(), 

53 station.strip(), 

54 location.strip(), 

55 channel.strip(), 

56 tmin, 

57 tmax, 

58 deltat, 

59 ydata) 

60 

61 tr.meta = { 

62 'offset_start': offset, 

63 'offset_end': tr_tuple[9], 

64 'last': tr_tuple[10], 

65 'segment_size': segment_size 

66 } 

67 

68 yield tr 

69 

70 if tr_tuple[10]: 

71 break 

72 

73 offset = tr_tuple[9] 

74 isegment += 1 

75 

76 except (OSError, mseed_ext.MSeedError) as e: 

77 raise FileLoadError(str(e)+' (file: %s)' % filename) 

78 

79 if have_zero_rate_traces: 

80 logger.warning( 

81 'Ignoring traces with sampling rate of zero in file %s ' 

82 '(maybe LOG traces)' % filename) 

83 

84 

85def as_tuple(tr, dataquality='D'): 

86 from pyrocko import mseed_ext 

87 itmin = int(round(tr.tmin*mseed_ext.HPTMODULUS)) 

88 itmax = int(round(tr.tmax*mseed_ext.HPTMODULUS)) 

89 srate = 1.0/tr.deltat 

90 return (tr.network, tr.station, tr.location, tr.channel, 

91 itmin, itmax, srate, dataquality, tr.get_ydata()) 

92 

93 

94def save(traces, filename_template, additional={}, overwrite=True, 

95 dataquality='D', record_length=4096, append=False, steim=1): 

96 from pyrocko import mseed_ext 

97 

98 assert record_length in VALID_RECORD_LENGTHS 

99 assert dataquality in ('D', 'E', 'C', 'O', 'T', 'L'), 'invalid dataquality' 

100 overwrite = True if append else overwrite 

101 

102 fn_tr = {} 

103 for tr in traces: 

104 for code, maxlen, val in zip( 

105 ['network', 'station', 'location', 'channel'], 

106 [2, 5, 2, 3], 

107 tr.nslc_id): 

108 

109 if len(val) > maxlen: 

110 raise CodeTooLong( 

111 '%s code too long to be stored in MSeed file: %s' % 

112 (code, val)) 

113 

114 fn = tr.fill_template(filename_template, **additional) 

115 if not overwrite and os.path.exists(fn): 

116 raise FileSaveError('File exists: %s' % fn) 

117 

118 if fn not in fn_tr: 

119 fn_tr[fn] = [] 

120 

121 fn_tr[fn].append(tr) 

122 

123 for fn, traces_thisfile in fn_tr.items(): 

124 trtups = [] 

125 traces_thisfile.sort(key=lambda a: a.full_id) 

126 for tr in traces_thisfile: 

127 trtups.append(as_tuple(tr, dataquality)) 

128 

129 ensuredirs(fn) 

130 try: 

131 mseed_ext.store_traces(trtups, fn, record_length, append, steim) 

132 except mseed_ext.MSeedError as e: 

133 raise FileSaveError( 

134 str(e) + " (while storing traces to file '%s')" % fn) 

135 

136 return list(fn_tr.keys()) 

137 

138 

139tcs = {} 

140 

141 

142def get_bytes(traces, dataquality='D', record_length=4096, steim=1): 

143 from pyrocko import mseed_ext 

144 

145 assert record_length in VALID_RECORD_LENGTHS 

146 assert dataquality in ('D', 'E', 'C', 'O', 'T', 'L'), 'invalid dataquality' 

147 

148 nbytes_approx = 0 

149 rl = record_length 

150 trtups = [] 

151 for tr in traces: 

152 for code, maxlen, val in zip( 

153 ['network', 'station', 'location', 'channel'], 

154 [2, 5, 2, 3], 

155 tr.nslc_id): 

156 

157 if len(val) > maxlen: 

158 raise CodeTooLong( 

159 '%s code too long to be stored in MSeed file: %s' % 

160 (code, val)) 

161 

162 nbytes_approx += math.ceil( 

163 tr.ydata.nbytes / (rl-MSEED_HEADER_BYTES)) * rl 

164 trtups.append(as_tuple(tr, dataquality)) 

165 

166 return mseed_ext.mseed_bytes(trtups, nbytes_approx, record_length, steim) 

167 

168 

169def detect(first512): 

170 

171 if len(first512) < 256: 

172 return False 

173 

174 rec = first512 

175 

176 try: 

177 sequence_number = int(rec[:6]) 

178 except Exception: 

179 return False 

180 if sequence_number < 0: 

181 return False 

182 

183 type_code = rec[6:7] 

184 if type_code in b'DRQM': 

185 bads = [] 

186 for sex in '<>': 

187 bad = False 

188 fmt = sex + '6s1s1s5s2s3s2s10sH2h4Bl2H' 

189 vals = unpack(fmt, rec[:48]) 

190 fmt_btime = sex + 'HHBBBBH' 

191 tvals = unpack(fmt_btime, vals[7]) 

192 if tvals[1] < 1 or tvals[1] > 367 or tvals[2] > 23 or \ 

193 tvals[3] > 59 or tvals[4] > 60 or tvals[6] > 9999: 

194 bad = True 

195 

196 bads.append(bad) 

197 

198 if all(bads): 

199 return False 

200 

201 else: 

202 if type_code not in b'VAST': 

203 return False 

204 

205 continuation_code = rec[7:8] 

206 if continuation_code != b' ': 

207 return False 

208 

209 blockette_type = rec[8:8+3].decode() 

210 if not re.match(r'^\d\d\d$', blockette_type): 

211 return False 

212 

213 try: 

214 blockette_length = int(rec[11:11+4]) 

215 except Exception: 

216 return False 

217 

218 if blockette_length < 7: 

219 return False 

220 

221 return True