1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6SAC IO library for Python
7'''
9import struct
10import logging
11import math
12import numpy as num
14from calendar import timegm
15from time import gmtime
17from pyrocko import trace, util
18from pyrocko.util import reuse
19from .io_common import FileLoadError
21logger = logging.getLogger('pyrocko.io.sac')
24def fixdoublefloat(x):
25 f = 10**math.floor(math.log10(x)) / 1000000.
26 return round(x/f)*f
29class SacError(Exception):
30 pass
33def nonetoempty(s):
34 if s is None:
35 return ''
36 else:
37 return s.strip()
40class SacFile(object):
41 nbytes_header = 632
42 header_num_format = {'little': '<70f40i', 'big': '>70f40i'}
44 header_keys = '''
45delta depmin depmax scale odelta b e o a internal0 t0 t1 t2 t3 t4 t5 t6 t7 t8
46t9 f resp0 resp1 resp2 resp3 resp4 resp5 resp6 resp7 resp8 resp9 stla stlo stel
47stdp evla evlo evel evdp mag user0 user1 user2 user3 user4 user5 user6 user7
48user8 user9 dist az baz gcarc internal1 internal2 depmen cmpaz cmpinc xminimum
49xmaximum yminimum ymaximum unused0 unused1 unused2 unused3 unused4 unused5
50unused6 nzyear nzjday nzhour nzmin nzsec nzmsec nvhdr norid nevid npts
51internal3 nwfid nxsize nysize unused7 iftype idep iztype unused8 iinst istreg
52ievreg ievtyp iqual isynth imagtyp imagsrc unused9 unused10 unused11 unused12
53unused13 unused14 unused15 unused16 leven lpspol lovrok lcalda unused17 kstnm
54kevnm khole ko ka kt0 kt1 kt2 kt3 kt4 kt5 kt6 kt7 kt8 kt9 kf kuser0 kuser1
55kuser2 kcmpnm knetwk kdatrd kinst
56'''.split()
58 header_enum_symbols = '''
59itime irlim iamph ixy iunkn idisp ivel iacc ib iday io ia it0 it1 it2 it3 it4
60it5 it6 it7 it8 it9 iradnv itannv iradev itanev inorth ieast ihorza idown iup
61illlbb iwwsn1 iwwsn2 ihglp isro inucl ipren ipostn iquake ipreq ipostq ichem
62iother igood iglch idrop ilowsn irldta ivolts ixyz imb ims iml imw imd imx
63ineic ipde iisc ireb iusgs ibrk icaltech illnl ievloc ijsop iuser iunknown iqb
64iqb1 iqb2 iqbx iqmt ieq ieq1 ieq2 ime iex inu inc io_ il ir it iu
65'''.split()
67 enum_header_vars = 'iftype idep iztype imagtype imagsrc ievtyp iqual ' \
68 'isynth'.split()
70 header_num2name = dict(
71 [(a+1, b) for (a, b) in enumerate(header_enum_symbols)])
72 header_name2num = dict(
73 [(b, a+1) for (a, b) in enumerate(header_enum_symbols)])
74 header_types = 'f'*70 + 'i'*35 + 'l'*5 + 'k'*23
75 undefined_value = {'f': -12345.0, 'i': -12345, 'l': None, 'k': '-12345'}
76 ldefaults = {
77 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'unused17': 0}
79 t_lookup = dict(zip(header_keys, header_types))
81 u_lookup = dict()
82 for k in header_keys:
83 u_lookup[k] = undefined_value[t_lookup[k]]
85 def ndatablocks(self):
86 '''
87 Get number of data blocks for this file's type.
88 '''
89 nblocks = {
90 'itime': 1, 'irlim': 2, 'iamph': 2, 'ixy': 2, 'ixyz': 3
91 }[SacFile.header_num2name[self.iftype]]
93 if nblocks == 1 and not self.leven:
94 nblocks = 2 # not sure about this...
96 return nblocks
98 def val_or_none(self, k, v):
99 '''
100 Replace SAC undef flags with None.
101 '''
102 if SacFile.u_lookup[k] == v:
103 return None
104 else:
105 return v
107 def get_ref_time(self):
108 '''
109 Get reference time as standard Unix timestamp.
110 '''
112 if None in (self.nzyear, self.nzjday, self.nzhour, self.nzmin,
113 self.nzsec, self.nzmsec):
114 raise SacError('Not all header values for reference time are set.')
116 return util.to_time_float(timegm(
117 (self.nzyear, 1, self.nzjday,
118 self.nzhour, self.nzmin, self.nzsec))) + self.nzmsec/1000.
120 def set_ref_time(self, timestamp):
121 '''
122 Set all header values defining reference time based on standard Unix
123 timestamp.
124 '''
126 secs = math.floor(timestamp)
127 msec = int(round((timestamp-secs)*1000.))
128 if msec == 1000:
129 secs += 1
130 msec = 0
132 t = gmtime(secs)
133 self.nzyear, self.nzjday, self.nzhour, self.nzmin, self.nzsec = \
134 t[0], t[7], t[3], t[4], t[5]
135 self.nzmsec = msec
137 def val_for_file(self, k, v):
138 '''
139 Convert attribute value to the form required when writing it to the
140 SAC file.
141 '''
143 t = SacFile.t_lookup[k]
144 if v is None:
145 if t == 'l':
146 return SacFile.ldefaults[k]
147 v = SacFile.u_lookup[k]
148 if t == 'f':
149 return float(v)
150 elif t == 'i':
151 return int(v)
152 elif t == 'l':
153 if v:
154 return 1
155 return 0
156 elif t == 'k':
157 ln = 8
158 if k == 'kevnm':
159 ln = 16 # only this header val has different length
160 return v.ljust(ln)[:ln]
162 def __init__(self, *args, **kwargs):
163 if 'from_trace' in kwargs:
164 self.clear()
165 trace = kwargs['from_trace']
166 if trace.meta:
167 for (k, v) in trace.meta.items():
168 if k in SacFile.header_keys:
169 setattr(self, k, v)
171 self.knetwk = trace.network
172 self.kstnm = trace.station
173 self.khole = trace.location
174 self.kcmpnm = trace.channel
175 self.set_ref_time(trace.tmin)
176 self.delta = trace.deltat
177 self.data = [trace.ydata.copy()]
178 self.npts = trace.ydata.size
179 self.b = 0.0
180 self.e = self.b + (self.npts-1)*self.delta
182 elif args:
183 self.read(*args, **kwargs)
184 else:
185 self.clear()
187 def clear(self):
188 '''
189 Empty file record.
190 '''
192 for k in SacFile.header_keys:
193 self.__dict__[k] = None
195 # set the required attributes
196 self.nvhdr = 6
197 self.iftype = SacFile.header_name2num['itime']
198 self.leven = True
199 self.delta = 1.0
200 self.npts = 0
201 self.b = 0.0
202 self.e = 0.0
203 self.data = [num.arange(0, dtype=num.float32)]
205 def check(self):
206 '''
207 Check the required header variables to have reasonable values.
208 '''
209 if self.iftype not in [SacFile.header_name2num[x] for x in
210 ('itime', 'irlim', 'iamph', 'ixy', 'ixyz')]:
211 raise SacError('Unknown SAC file type: %i.' % self.iftype)
212 if self.nvhdr < 1 or 20 < self.nvhdr:
213 raise SacError('Unreasonable SAC header version number found.')
214 if self.npts < 0:
215 raise SacError(
216 'Negative number of samples specified in NPTS header.')
217 if self.leven not in (0, 1, -12345):
218 raise SacError('Header value LEVEN must be either 0 or 1.')
219 if self.leven and self.delta <= 0.0:
220 raise SacError(
221 'Header value DELTA should be positive for evenly spaced '
222 'samples')
223 if self.e is not None and self.b > self.e:
224 raise SacError(
225 'Beginning value of independent variable greater than its '
226 'ending value.')
227 if self.nvhdr != 6:
228 logging.warn(
229 'This module has only been tested with SAC header version 6.'
230 'This file has header version %i. '
231 'It might still work though...' % self.nvhdr)
233 def read(self, filename, load_data=True, byte_sex='try'):
234 '''
235 Read SAC file.
237 filename -- Name of SAC file.
238 load_data -- If True, the data is read, otherwise only read headers.
239 byte_sex -- Endianness: 'try', 'little' or 'big'
240 '''
242 nbh = SacFile.nbytes_header
244 # read in all data
245 with open(filename, 'rb') as f:
246 if load_data:
247 filedata = f.read()
248 else:
249 filedata = f.read(nbh)
251 if len(filedata) < nbh:
252 raise SacError('File too short to be a SAC file: %s' % filename)
254 # possibly try out different endiannesses
255 if byte_sex == 'try':
256 sexes = ('little', 'big')
257 else:
258 sexes = (byte_sex,)
260 for isex, sex in enumerate(sexes):
261 format = SacFile.header_num_format[sex]
262 nbn = struct.calcsize(format)
263 hv = list(struct.unpack(format, filedata[:nbn]))
265 strings = str(filedata[nbn:nbh].decode('ascii'))
266 hv.append(strings[:8].rstrip(' \x00'))
267 hv.append(strings[8:24].rstrip(' \x00'))
268 for i in range(len(strings[24:])//8):
269 hv.append(strings[24+i*8:24+(i+1)*8].rstrip(' \x00'))
271 self.header_vals = hv
272 for k, v in zip(SacFile.header_keys, self.header_vals):
273 vn = self.val_or_none(k, v)
274 self.__dict__[k] = vn
276 if self.leven == -12345:
277 self.leven = True
279 self.data = []
280 try:
281 self.check()
282 break
283 except SacError as e:
284 if isex == len(sexes)-1:
285 raise e
287 self.delta = fixdoublefloat(self.delta)
289 if byte_sex == 'try':
290 logger.debug(
291 'This seems to be a %s endian SAC file: %s' % (sex, filename))
293 # possibly get data
294 if load_data:
295 nblocks = self.ndatablocks()
296 nbb = self.npts*4 # word length is always 4 bytes in sac files
297 for iblock in range(nblocks):
298 if len(filedata) < nbh+(iblock+1)*nbb:
299 raise SacError('File is incomplete.')
301 if sex == 'big':
302 dtype = num.dtype('>f4')
303 else:
304 dtype = num.dtype('<f4')
306 self.data.append(num.array(num.frombuffer(
307 filedata[nbh+iblock*nbb:nbh+(iblock+1)*nbb],
308 dtype=dtype),
309 dtype=float))
311 if len(filedata) > nbh+nblocks*nbb:
312 logger.warning(
313 'Unused data (%i bytes) at end of SAC file: %s (npts=%i)'
314 % (len(filedata) - nbh+nblocks*nbb, filename, self.npts))
316 def write(self, filename, byte_sex='little'):
317 '''
318 Write SAC file.
319 '''
321 self.check()
323 # create header data
324 format = SacFile.header_num_format[byte_sex]
325 numerical_values = []
326 string_values = []
327 for k in SacFile.header_keys:
328 v = self.__dict__[k]
329 vv = self.val_for_file(k, v)
330 if SacFile.t_lookup[k] == 'k':
331 string_values.append(vv)
332 else:
333 numerical_values.append(vv)
335 header_data = struct.pack(format, *numerical_values)
336 header_data += bytes(''.join(string_values).encode('ascii'))
338 # check that enough data is available
339 nblocks = self.ndatablocks()
340 if len(self.data) != nblocks:
341 raise SacError(
342 'Need %i data blocks for file type %s.'
343 % (nblocks, SacFile.header_num2name[self.iftype]))
345 for fdata in self.data:
346 if len(fdata) != self.npts:
347 raise SacError(
348 'Data length (%i) does not match NPTS header value (%i)'
349 % (len(fdata), self.npts))
351 # dump data to file
352 with open(filename, 'wb') as f:
353 f.write(header_data)
354 for fdata in self.data:
355 f.write(fdata.astype(num.float32).tobytes())
357 def __str__(self):
358 str = ''
359 for k in SacFile.header_keys:
360 v = self.__dict__[k]
361 if v is not None:
362 if k in SacFile.enum_header_vars:
363 if v in SacFile.header_num2name:
364 v = SacFile.header_num2name[v]
365 str += '%s: %s\n' % (k, v)
367 return str
369 def to_trace(self):
371 assert self.iftype == SacFile.header_name2num['itime']
372 assert self.leven
374 tmin = self.get_ref_time() + self.b
375 tmax = tmin + self.delta*(self.npts-1)
377 data = None
378 if self.data:
379 data = self.data[0]
381 meta = {}
382 exclude = ('b', 'e', 'knetwk', 'kstnm', 'khole', 'kcmpnm', 'delta',
383 'nzyear', 'nzjday', 'nzhour', 'nzmin', 'nzsec', 'nzmsec')
385 for k in SacFile.header_keys:
386 if k in exclude:
387 continue
388 v = self.__dict__[k]
389 if v is not None:
390 meta[reuse(k)] = v
392 return trace.Trace(
393 nonetoempty(self.knetwk)[:2],
394 nonetoempty(self.kstnm)[:5],
395 nonetoempty(self.khole)[:2],
396 nonetoempty(self.kcmpnm)[:3],
397 tmin,
398 tmax,
399 self.delta,
400 data,
401 meta=meta)
404def iload(filename, load_data=True):
406 try:
407 sacf = SacFile(filename, load_data=load_data)
408 tr = sacf.to_trace()
409 yield tr
411 except (OSError, SacError) as e:
412 raise FileLoadError(e)
415def detect(first512):
417 if len(first512) < 512: # SAC header is 632 bytes long
418 return False
420 for sex in 'little', 'big':
421 format = SacFile.header_num_format[sex]
422 nbn = struct.calcsize(format)
424 hv = list(struct.unpack(format, first512[:nbn]))
425 iftype, nvhdr, npts, leven, delta, e, b = [
426 hv[i] for i in (85, 76, 79, 105, 0, 6, 5)]
428 if iftype not in [SacFile.header_name2num[x] for x in (
429 'itime', 'irlim', 'iamph', 'ixy', 'ixyz')]:
430 continue
431 if nvhdr < 1 or 20 < nvhdr:
432 continue
433 if npts < 0:
434 continue
435 if leven not in (0, 1, -12345):
436 continue
437 if leven and delta <= 0.0:
438 continue
439 if e != -12345.0 and b > e:
440 continue
442 return True
444 return False