1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6SAC IO library for Python 

7''' 

8from __future__ import absolute_import 

9 

10import struct 

11import logging 

12import math 

13import numpy as num 

14 

15from calendar import timegm 

16from time import gmtime 

17 

18from pyrocko import trace, util 

19from pyrocko.util import reuse 

20from .io_common import FileLoadError 

21 

22logger = logging.getLogger('pyrocko.io.sac') 

23 

24 

25def fixdoublefloat(x): 

26 f = 10**math.floor(math.log10(x)) / 1000000. 

27 return round(x/f)*f 

28 

29 

30class SacError(Exception): 

31 pass 

32 

33 

34def nonetoempty(s): 

35 if s is None: 

36 return '' 

37 else: 

38 return s.strip() 

39 

40 

41class SacFile(object): 

42 nbytes_header = 632 

43 header_num_format = {'little': '<70f40i', 'big': '>70f40i'} 

44 

45 header_keys = ''' 

46delta depmin depmax scale odelta b e o a internal0 t0 t1 t2 t3 t4 t5 t6 t7 t8 

47t9 f resp0 resp1 resp2 resp3 resp4 resp5 resp6 resp7 resp8 resp9 stla stlo stel 

48stdp evla evlo evel evdp mag user0 user1 user2 user3 user4 user5 user6 user7 

49user8 user9 dist az baz gcarc internal1 internal2 depmen cmpaz cmpinc xminimum 

50xmaximum yminimum ymaximum unused0 unused1 unused2 unused3 unused4 unused5 

51unused6 nzyear nzjday nzhour nzmin nzsec nzmsec nvhdr norid nevid npts 

52internal3 nwfid nxsize nysize unused7 iftype idep iztype unused8 iinst istreg 

53ievreg ievtyp iqual isynth imagtyp imagsrc unused9 unused10 unused11 unused12 

54unused13 unused14 unused15 unused16 leven lpspol lovrok lcalda unused17 kstnm 

55kevnm khole ko ka kt0 kt1 kt2 kt3 kt4 kt5 kt6 kt7 kt8 kt9 kf kuser0 kuser1 

56kuser2 kcmpnm knetwk kdatrd kinst 

57'''.split() 

58 

59 header_enum_symbols = ''' 

60itime irlim iamph ixy iunkn idisp ivel iacc ib iday io ia it0 it1 it2 it3 it4 

61it5 it6 it7 it8 it9 iradnv itannv iradev itanev inorth ieast ihorza idown iup 

62illlbb iwwsn1 iwwsn2 ihglp isro inucl ipren ipostn iquake ipreq ipostq ichem 

63iother igood iglch idrop ilowsn irldta ivolts ixyz imb ims iml imw imd imx 

64ineic ipde iisc ireb iusgs ibrk icaltech illnl ievloc ijsop iuser iunknown iqb 

65iqb1 iqb2 iqbx iqmt ieq ieq1 ieq2 ime iex inu inc io_ il ir it iu 

66'''.split() 

67 

68 enum_header_vars = 'iftype idep iztype imagtype imagsrc ievtyp iqual ' \ 

69 'isynth'.split() 

70 

71 header_num2name = dict( 

72 [(a+1, b) for (a, b) in enumerate(header_enum_symbols)]) 

73 header_name2num = dict( 

74 [(b, a+1) for (a, b) in enumerate(header_enum_symbols)]) 

75 header_types = 'f'*70 + 'i'*35 + 'l'*5 + 'k'*23 

76 undefined_value = {'f': -12345.0, 'i': -12345, 'l': None, 'k': '-12345'} 

77 ldefaults = { 

78 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'unused17': 0} 

79 

80 t_lookup = dict(zip(header_keys, header_types)) 

81 

82 u_lookup = dict() 

83 for k in header_keys: 

84 u_lookup[k] = undefined_value[t_lookup[k]] 

85 

86 def ndatablocks(self): 

87 ''' 

88 Get number of data blocks for this file's type. 

89 ''' 

90 nblocks = { 

91 'itime': 1, 'irlim': 2, 'iamph': 2, 'ixy': 2, 'ixyz': 3 

92 }[SacFile.header_num2name[self.iftype]] 

93 

94 if nblocks == 1 and not self.leven: 

95 nblocks = 2 # not sure about this... 

96 

97 return nblocks 

98 

99 def val_or_none(self, k, v): 

100 ''' 

101 Replace SAC undef flags with None. 

102 ''' 

103 if SacFile.u_lookup[k] == v: 

104 return None 

105 else: 

106 return v 

107 

108 def get_ref_time(self): 

109 ''' 

110 Get reference time as standard Unix timestamp. 

111 ''' 

112 

113 if None in (self.nzyear, self.nzjday, self.nzhour, self.nzmin, 

114 self.nzsec, self.nzmsec): 

115 raise SacError('Not all header values for reference time are set.') 

116 

117 return util.to_time_float(timegm( 

118 (self.nzyear, 1, self.nzjday, 

119 self.nzhour, self.nzmin, self.nzsec))) + self.nzmsec/1000. 

120 

121 def set_ref_time(self, timestamp): 

122 ''' 

123 Set all header values defining reference time based on standard Unix 

124 timestamp. 

125 ''' 

126 

127 secs = math.floor(timestamp) 

128 msec = int(round((timestamp-secs)*1000.)) 

129 if msec == 1000: 

130 secs += 1 

131 msec = 0 

132 

133 t = gmtime(secs) 

134 self.nzyear, self.nzjday, self.nzhour, self.nzmin, self.nzsec = \ 

135 t[0], t[7], t[3], t[4], t[5] 

136 self.nzmsec = msec 

137 

138 def val_for_file(self, k, v): 

139 ''' 

140 Convert attribute value to the form required when writing it to the 

141 SAC file. 

142 ''' 

143 

144 t = SacFile.t_lookup[k] 

145 if v is None: 

146 if t == 'l': 

147 return SacFile.ldefaults[k] 

148 v = SacFile.u_lookup[k] 

149 if t == 'f': 

150 return float(v) 

151 elif t == 'i': 

152 return int(v) 

153 elif t == 'l': 

154 if v: 

155 return 1 

156 return 0 

157 elif t == 'k': 

158 ln = 8 

159 if k == 'kevnm': 

160 ln = 16 # only this header val has different length 

161 return v.ljust(ln)[:ln] 

162 

163 def __init__(self, *args, **kwargs): 

164 if 'from_trace' in kwargs: 

165 self.clear() 

166 trace = kwargs['from_trace'] 

167 if trace.meta: 

168 for (k, v) in trace.meta.items(): 

169 if k in SacFile.header_keys: 

170 setattr(self, k, v) 

171 

172 self.knetwk = trace.network 

173 self.kstnm = trace.station 

174 self.khole = trace.location 

175 self.kcmpnm = trace.channel 

176 self.set_ref_time(trace.tmin) 

177 self.delta = trace.deltat 

178 self.data = [trace.ydata.copy()] 

179 self.npts = trace.ydata.size 

180 self.b = 0.0 

181 self.e = self.b + (self.npts-1)*self.delta 

182 

183 elif args: 

184 self.read(*args, **kwargs) 

185 else: 

186 self.clear() 

187 

188 def clear(self): 

189 ''' 

190 Empty file record. 

191 ''' 

192 

193 for k in SacFile.header_keys: 

194 self.__dict__[k] = None 

195 

196 # set the required attributes 

197 self.nvhdr = 6 

198 self.iftype = SacFile.header_name2num['itime'] 

199 self.leven = True 

200 self.delta = 1.0 

201 self.npts = 0 

202 self.b = 0.0 

203 self.e = 0.0 

204 self.data = [num.arange(0, dtype=num.float32)] 

205 

206 def check(self): 

207 ''' 

208 Check the required header variables to have reasonable values. 

209 ''' 

210 if self.iftype not in [SacFile.header_name2num[x] for x in 

211 ('itime', 'irlim', 'iamph', 'ixy', 'ixyz')]: 

212 raise SacError('Unknown SAC file type: %i.' % self.iftype) 

213 if self.nvhdr < 1 or 20 < self.nvhdr: 

214 raise SacError('Unreasonable SAC header version number found.') 

215 if self.npts < 0: 

216 raise SacError( 

217 'Negative number of samples specified in NPTS header.') 

218 if self.leven not in (0, 1, -12345): 

219 raise SacError('Header value LEVEN must be either 0 or 1.') 

220 if self.leven and self.delta <= 0.0: 

221 raise SacError( 

222 'Header value DELTA should be positive for evenly spaced ' 

223 'samples') 

224 if self.e is not None and self.b > self.e: 

225 raise SacError( 

226 'Beginning value of independent variable greater than its ' 

227 'ending value.') 

228 if self.nvhdr != 6: 

229 logging.warn( 

230 'This module has only been tested with SAC header version 6.' 

231 'This file has header version %i. ' 

232 'It might still work though...' % self.nvhdr) 

233 

234 def read(self, filename, load_data=True, byte_sex='try'): 

235 ''' 

236 Read SAC file. 

237 

238 filename -- Name of SAC file. 

239 load_data -- If True, the data is read, otherwise only read headers. 

240 byte_sex -- Endianness: 'try', 'little' or 'big' 

241 ''' 

242 

243 nbh = SacFile.nbytes_header 

244 

245 # read in all data 

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

247 if load_data: 

248 filedata = f.read() 

249 else: 

250 filedata = f.read(nbh) 

251 

252 if len(filedata) < nbh: 

253 raise SacError('File too short to be a SAC file: %s' % filename) 

254 

255 # possibly try out different endiannesses 

256 if byte_sex == 'try': 

257 sexes = ('little', 'big') 

258 else: 

259 sexes = (byte_sex,) 

260 

261 for isex, sex in enumerate(sexes): 

262 format = SacFile.header_num_format[sex] 

263 nbn = struct.calcsize(format) 

264 hv = list(struct.unpack(format, filedata[:nbn])) 

265 

266 strings = str(filedata[nbn:nbh].decode('ascii')) 

267 hv.append(strings[:8].rstrip(' \x00')) 

268 hv.append(strings[8:24].rstrip(' \x00')) 

269 for i in range(len(strings[24:])//8): 

270 hv.append(strings[24+i*8:24+(i+1)*8].rstrip(' \x00')) 

271 

272 self.header_vals = hv 

273 for k, v in zip(SacFile.header_keys, self.header_vals): 

274 vn = self.val_or_none(k, v) 

275 self.__dict__[k] = vn 

276 

277 if self.leven == -12345: 

278 self.leven = True 

279 

280 self.data = [] 

281 try: 

282 self.check() 

283 break 

284 except SacError as e: 

285 if isex == len(sexes)-1: 

286 raise e 

287 

288 self.delta = fixdoublefloat(self.delta) 

289 

290 if byte_sex == 'try': 

291 logger.debug( 

292 'This seems to be a %s endian SAC file: %s' % (sex, filename)) 

293 

294 # possibly get data 

295 if load_data: 

296 nblocks = self.ndatablocks() 

297 nbb = self.npts*4 # word length is always 4 bytes in sac files 

298 for iblock in range(nblocks): 

299 if len(filedata) < nbh+(iblock+1)*nbb: 

300 raise SacError('File is incomplete.') 

301 

302 if sex == 'big': 

303 dtype = num.dtype('>f4') 

304 else: 

305 dtype = num.dtype('<f4') 

306 

307 self.data.append(num.array(num.frombuffer( 

308 filedata[nbh+iblock*nbb:nbh+(iblock+1)*nbb], 

309 dtype=dtype), 

310 dtype=float)) 

311 

312 if len(filedata) > nbh+nblocks*nbb: 

313 logger.warning( 

314 'Unused data (%i bytes) at end of SAC file: %s (npts=%i)' 

315 % (len(filedata) - nbh+nblocks*nbb, filename, self.npts)) 

316 

317 def write(self, filename, byte_sex='little'): 

318 ''' 

319 Write SAC file. 

320 ''' 

321 

322 self.check() 

323 

324 # create header data 

325 format = SacFile.header_num_format[byte_sex] 

326 numerical_values = [] 

327 string_values = [] 

328 for k in SacFile.header_keys: 

329 v = self.__dict__[k] 

330 vv = self.val_for_file(k, v) 

331 if SacFile.t_lookup[k] == 'k': 

332 string_values.append(vv) 

333 else: 

334 numerical_values.append(vv) 

335 

336 header_data = struct.pack(format, *numerical_values) 

337 header_data += bytes(''.join(string_values).encode('ascii')) 

338 

339 # check that enough data is available 

340 nblocks = self.ndatablocks() 

341 if len(self.data) != nblocks: 

342 raise SacError( 

343 'Need %i data blocks for file type %s.' 

344 % (nblocks, SacFile.header_num2name[self.iftype])) 

345 

346 for fdata in self.data: 

347 if len(fdata) != self.npts: 

348 raise SacError( 

349 'Data length (%i) does not match NPTS header value (%i)' 

350 % (len(fdata), self.npts)) 

351 

352 # dump data to file 

353 with open(filename, 'wb') as f: 

354 f.write(header_data) 

355 for fdata in self.data: 

356 f.write(fdata.astype(num.float32).tobytes()) 

357 

358 def __str__(self): 

359 str = '' 

360 for k in SacFile.header_keys: 

361 v = self.__dict__[k] 

362 if v is not None: 

363 if k in SacFile.enum_header_vars: 

364 if v in SacFile.header_num2name: 

365 v = SacFile.header_num2name[v] 

366 str += '%s: %s\n' % (k, v) 

367 

368 return str 

369 

370 def to_trace(self): 

371 

372 assert self.iftype == SacFile.header_name2num['itime'] 

373 assert self.leven 

374 

375 tmin = self.get_ref_time() + self.b 

376 tmax = tmin + self.delta*(self.npts-1) 

377 

378 data = None 

379 if self.data: 

380 data = self.data[0] 

381 

382 meta = {} 

383 exclude = ('b', 'e', 'knetwk', 'kstnm', 'khole', 'kcmpnm', 'delta', 

384 'nzyear', 'nzjday', 'nzhour', 'nzmin', 'nzsec', 'nzmsec') 

385 

386 for k in SacFile.header_keys: 

387 if k in exclude: 

388 continue 

389 v = self.__dict__[k] 

390 if v is not None: 

391 meta[reuse(k)] = v 

392 

393 return trace.Trace( 

394 nonetoempty(self.knetwk)[:2], 

395 nonetoempty(self.kstnm)[:5], 

396 nonetoempty(self.khole)[:2], 

397 nonetoempty(self.kcmpnm)[:3], 

398 tmin, 

399 tmax, 

400 self.delta, 

401 data, 

402 meta=meta) 

403 

404 

405def iload(filename, load_data=True): 

406 

407 try: 

408 sacf = SacFile(filename, load_data=load_data) 

409 tr = sacf.to_trace() 

410 yield tr 

411 

412 except (OSError, SacError) as e: 

413 raise FileLoadError(e) 

414 

415 

416def detect(first512): 

417 

418 if len(first512) < 512: # SAC header is 632 bytes long 

419 return False 

420 

421 for sex in 'little', 'big': 

422 format = SacFile.header_num_format[sex] 

423 nbn = struct.calcsize(format) 

424 

425 hv = list(struct.unpack(format, first512[:nbn])) 

426 iftype, nvhdr, npts, leven, delta, e, b = [ 

427 hv[i] for i in (85, 76, 79, 105, 0, 6, 5)] 

428 

429 if iftype not in [SacFile.header_name2num[x] for x in ( 

430 'itime', 'irlim', 'iamph', 'ixy', 'ixyz')]: 

431 continue 

432 if nvhdr < 1 or 20 < nvhdr: 

433 continue 

434 if npts < 0: 

435 continue 

436 if leven not in (0, 1, -12345): 

437 continue 

438 if leven and delta <= 0.0: 

439 continue 

440 if e != -12345.0 and b > e: 

441 continue 

442 

443 return True 

444 

445 return False