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

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Reader for the `Güralp GCF 

8<https://www.guralp.com/customer-support/common-questions/other-data-formats/data-formats/gcf>`_ 

9format. 

10''' 

11 

12import sys 

13import struct 

14import re 

15import numpy as num 

16 

17from io import StringIO 

18from collections import namedtuple 

19 

20from pyrocko import util, trace 

21from .io_common import FileLoadError 

22 

23 

24g_guralp_zero = None 

25 

26 

27def get_guralp_zero(): 

28 global g_guralp_zero 

29 if g_guralp_zero is None: 

30 g_guralp_zero = util.str_to_time('1989-11-17 00:00:00') 

31 

32 return g_guralp_zero 

33 

34 

35class GCFLoadError(FileLoadError): 

36 pass 

37 

38 

39class EOF(Exception): 

40 pass 

41 

42 

43Header = namedtuple( 

44 'Header', 

45 ''' 

46 block_type 

47 system_id 

48 stream_id 

49 instrument_type 

50 time 

51 gain 

52 ttl 

53 sample_rate 

54 compression 

55 nrecords 

56 '''.split()) 

57 

58 

59def read(f, n, eof_ok=False): 

60 data = f.read(n) 

61 if eof_ok and len(data) == 0: 

62 raise EOF() 

63 

64 if len(data) != n: 

65 raise GCFLoadError('Unexpected end of file') 

66 

67 return data 

68 

69 

70def read_header(f, endianness='>'): 

71 e = endianness 

72 data = read(f, 16, eof_ok=True) 

73 

74 isystem_id, istream_id = struct.unpack(e+'II', data[:8]) 

75 ex = isystem_id & 0x80000000 

76 if ex: 

77 ex2 = isystem_id & (1 << 30) 

78 if ex2: 

79 system_id = util.base36encode(isystem_id & (2**21 - 1)) 

80 else: 

81 system_id = util.base36encode(isystem_id & (2**26 - 1)) 

82 

83 instrument_type = (isystem_id >> 26) & 0b1 

84 gain = [None, 1, 2, 4, 8, 16, 32, 64][(isystem_id >> 27) & 0b111] 

85 else: 

86 

87 system_id = util.base36encode(isystem_id) 

88 instrument_type = None 

89 gain = None 

90 

91 stream_id = util.base36encode(istream_id) 

92 

93 i_day_second = struct.unpack(e+'I', data[8:12])[0] 

94 iday = i_day_second >> 17 

95 isecond = i_day_second & 0x1ffff 

96 time = (iday*24*60*60) + get_guralp_zero() + isecond 

97 

98 ittl, israte, compression, nrecords = struct.unpack(e+'BBBB', data[12:]) 

99 if nrecords > 250: 

100 raise FileLoadError('Header indicates too many records in block.') 

101 

102 if israte == 0: 

103 if compression == 4 and stream_id[-2:] == '00': 

104 block_type = 'status_block' 

105 

106 elif compression == 4 and stream_id[-2:] == '01': 

107 block_type = 'unified_status_block' 

108 

109 elif compression == 4 and stream_id[-2:] == 'SM': 

110 block_type = 'strong_motion_block' 

111 

112 elif stream_id[-2:] == 'CD': 

113 block_type = 'cd_status_block' 

114 

115 elif compression == 4 and stream_id[-2:] == 'BP': 

116 block_type = 'byte_pipe_block' 

117 

118 elif stream_id[-2:] == 'IB': 

119 block_type = 'information_block' 

120 

121 else: 

122 raise FileLoadError('Unexpected block type found.') 

123 

124 return Header(block_type, system_id, stream_id, instrument_type, 

125 time, gain, ittl, 0.0, compression, nrecords) 

126 else: 

127 block_type = 'data_block' 

128 

129 if not re.match(r'^([ZNEXC][0-9A-CG-S]|M[0-9A-F])$', stream_id[-2:]): 

130 raise FileLoadError('Unexpected data stream ID') 

131 

132 sample_rate_tab = { 

133 157: (0.1, None), 

134 161: (0.125, None), 

135 162: (0.2, None), 

136 164: (0.25, None), 

137 167: (0.5, None), 

138 171: (400., 8), 

139 174: (500., 2), 

140 176: (1000., 4), 

141 179: (2000., 8), 

142 181: (4000., 16)} 

143 

144 if israte in sample_rate_tab: 

145 sample_rate, tfod = sample_rate_tab[israte] 

146 else: 

147 sample_rate = float(israte) 

148 tfod = None 

149 

150 if tfod is not None: 

151 toff = (compression >> 4) // tfod 

152 compression = compression & 0b1111 

153 time += toff 

154 

155 if compression not in (1, 2, 4): 

156 raise GCFLoadError( 

157 'Unsupported compression code: %i' % compression) 

158 

159 return Header(block_type, system_id, stream_id, instrument_type, 

160 time, gain, ittl, sample_rate, compression, nrecords) 

161 

162 

163def read_data(f, h, endianness='>'): 

164 e = endianness 

165 data = read(f, 1024 - 16) 

166 first = struct.unpack(e+'i', data[0:4])[0] 

167 dtype = {1: e+'i4', 2: e+'i2', 4: e+'i1'} 

168 if h.compression not in dtype: 

169 raise GCFLoadError('Unsupported compression code: %i' % h.compression) 

170 

171 nsamples = h.compression * h.nrecords 

172 difs = num.frombuffer(data[4:4+h.nrecords*4], dtype=dtype[h.compression], 

173 count=nsamples) 

174 samples = difs.astype(num.int32) 

175 samples[0] += first 

176 samples = num.cumsum(samples, dtype=num.int32) 

177 last = struct.unpack(e+'i', data[4+h.nrecords*4:4+h.nrecords*4+4])[0] 

178 if last != samples[-1]: 

179 raise GCFLoadError('Checksum error occured') 

180 

181 return samples 

182 

183 

184def read_status(f, h): 

185 data = read(f, 1024 - 16) 

186 return data[:h.nrecords*4] 

187 

188 

189def iload(filename, load_data=True): 

190 traces = {} 

191 

192 with open(filename, 'rb') as f: 

193 try: 

194 while True: 

195 h = read_header(f) 

196 if h.block_type == 'data_block': 

197 deltat = 1.0 / h.sample_rate 

198 if load_data: 

199 samples = read_data(f, h) 

200 tmax = None 

201 else: 

202 f.seek(1024 - 16, 1) 

203 samples = None 

204 tmax = h.time + ( 

205 h.nrecords * h.compression - 1) * deltat 

206 

207 nslc = ('', h.system_id, '', h.stream_id) 

208 

209 if nslc in traces: 

210 tr = traces[nslc] 

211 if abs((tr.tmax + tr.deltat) - h.time) < deltat*0.0001: 

212 if samples is not None: 

213 tr.append(samples) 

214 else: 

215 tr.tmax = tmax 

216 

217 else: 

218 del traces[nslc] 

219 yield tr 

220 

221 if nslc not in traces: 

222 traces[nslc] = trace.Trace( 

223 *nslc, 

224 tmin=h.time, 

225 deltat=deltat, 

226 ydata=samples, 

227 tmax=tmax) 

228 

229 else: 

230 f.seek(1024 - 16, 1) 

231 

232 except EOF: 

233 for tr in traces.values(): 

234 yield tr 

235 

236 

237def detect(first512): 

238 # does not work properly, produces too many false positives 

239 # difficult to improve due to the very compact GCF header 

240 try: 

241 if len(first512) < 512: 

242 return False 

243 

244 f = StringIO(first512) 

245 read_header(f) 

246 return True 

247 

248 except Exception: 

249 return False 

250 

251 

252if __name__ == '__main__': 

253 util.setup_logging('warn') 

254 

255 all_traces = [] 

256 for fn in sys.argv[1:]: 

257 if detect(open(fn, 'rb').read(512)): 

258 print(fn) 

259 all_traces.extend(iload(fn)) 

260 

261 trace.snuffle(all_traces)