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
10import struct
11import logging
12import math
13import numpy as num
15from calendar import timegm
16from time import gmtime
18from pyrocko import trace, util
19from pyrocko.util import reuse
20from .io_common import FileLoadError
22logger = logging.getLogger('pyrocko.io.sac')
25def fixdoublefloat(x):
26 f = 10**math.floor(math.log10(x)) / 1000000.
27 return round(x/f)*f
30class SacError(Exception):
31 pass
34def nonetoempty(s):
35 if s is None:
36 return ''
37 else:
38 return s.strip()
41class SacFile(object):
42 nbytes_header = 632
43 header_num_format = {'little': '<70f40i', 'big': '>70f40i'}
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()
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()
68 enum_header_vars = 'iftype idep iztype imagtype imagsrc ievtyp iqual ' \
69 'isynth'.split()
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}
80 t_lookup = dict(zip(header_keys, header_types))
82 u_lookup = dict()
83 for k in header_keys:
84 u_lookup[k] = undefined_value[t_lookup[k]]
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]]
94 if nblocks == 1 and not self.leven:
95 nblocks = 2 # not sure about this...
97 return nblocks
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
108 def get_ref_time(self):
109 '''
110 Get reference time as standard Unix timestamp.
111 '''
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.')
117 return util.to_time_float(timegm(
118 (self.nzyear, 1, self.nzjday,
119 self.nzhour, self.nzmin, self.nzsec))) + self.nzmsec/1000.
121 def set_ref_time(self, timestamp):
122 '''
123 Set all header values defining reference time based on standard Unix
124 timestamp.
125 '''
127 secs = math.floor(timestamp)
128 msec = int(round((timestamp-secs)*1000.))
129 if msec == 1000:
130 secs += 1
131 msec = 0
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
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 '''
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]
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)
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
183 elif args:
184 self.read(*args, **kwargs)
185 else:
186 self.clear()
188 def clear(self):
189 '''
190 Empty file record.
191 '''
193 for k in SacFile.header_keys:
194 self.__dict__[k] = None
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)]
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)
234 def read(self, filename, load_data=True, byte_sex='try'):
235 '''
236 Read SAC file.
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 '''
243 nbh = SacFile.nbytes_header
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)
252 if len(filedata) < nbh:
253 raise SacError('File too short to be a SAC file: %s' % filename)
255 # possibly try out different endiannesses
256 if byte_sex == 'try':
257 sexes = ('little', 'big')
258 else:
259 sexes = (byte_sex,)
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]))
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'))
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
277 if self.leven == -12345:
278 self.leven = True
280 self.data = []
281 try:
282 self.check()
283 break
284 except SacError as e:
285 if isex == len(sexes)-1:
286 raise e
288 self.delta = fixdoublefloat(self.delta)
290 if byte_sex == 'try':
291 logger.debug(
292 'This seems to be a %s endian SAC file: %s' % (sex, filename))
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.')
302 if sex == 'big':
303 dtype = num.dtype('>f4')
304 else:
305 dtype = num.dtype('<f4')
307 self.data.append(num.array(num.frombuffer(
308 filedata[nbh+iblock*nbb:nbh+(iblock+1)*nbb],
309 dtype=dtype),
310 dtype=float))
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))
317 def write(self, filename, byte_sex='little'):
318 '''
319 Write SAC file.
320 '''
322 self.check()
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)
336 header_data = struct.pack(format, *numerical_values)
337 header_data += bytes(''.join(string_values).encode('ascii'))
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]))
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))
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())
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)
368 return str
370 def to_trace(self):
372 assert self.iftype == SacFile.header_name2num['itime']
373 assert self.leven
375 tmin = self.get_ref_time() + self.b
376 tmax = tmin + self.delta*(self.npts-1)
378 data = None
379 if self.data:
380 data = self.data[0]
382 meta = {}
383 exclude = ('b', 'e', 'knetwk', 'kstnm', 'khole', 'kcmpnm', 'delta',
384 'nzyear', 'nzjday', 'nzhour', 'nzmin', 'nzsec', 'nzmsec')
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
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)
405def iload(filename, load_data=True):
407 try:
408 sacf = SacFile(filename, load_data=load_data)
409 tr = sacf.to_trace()
410 yield tr
412 except (OSError, SacError) as e:
413 raise FileLoadError(e)
416def detect(first512):
418 if len(first512) < 512: # SAC header is 632 bytes long
419 return False
421 for sex in 'little', 'big':
422 format = SacFile.header_num_format[sex]
423 nbn = struct.calcsize(format)
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)]
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
443 return True
445 return False