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())
268 hv.append(strings[8:24].rstrip())
269 for i in range(len(strings[24:])//8):
270 hv.append(strings[24+i*8:24+(i+1)*8].rstrip())
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 if isinstance(vn, str):
276 ipos = vn.find('\x00')
277 if ipos != -1:
278 vn = vn[:ipos]
280 self.__dict__[k] = vn.rstrip()
281 else:
282 self.__dict__[k] = vn
284 if self.leven == -12345:
285 self.leven = True
287 self.data = []
288 try:
289 self.check()
290 break
291 except SacError as e:
292 if isex == len(sexes)-1:
293 raise e
295 self.delta = fixdoublefloat(self.delta)
297 if byte_sex == 'try':
298 logger.debug(
299 'This seems to be a %s endian SAC file: %s' % (sex, filename))
301 # possibly get data
302 if load_data:
303 nblocks = self.ndatablocks()
304 nbb = self.npts*4 # word length is always 4 bytes in sac files
305 for iblock in range(nblocks):
306 if len(filedata) < nbh+(iblock+1)*nbb:
307 raise SacError('File is incomplete.')
309 if sex == 'big':
310 dtype = num.dtype('>f4')
311 else:
312 dtype = num.dtype('<f4')
314 self.data.append(num.array(num.fromstring(
315 filedata[nbh+iblock*nbb:nbh+(iblock+1)*nbb],
316 dtype=dtype),
317 dtype=float))
319 if len(filedata) > nbh+nblocks*nbb:
320 logger.warning(
321 'Unused data (%i bytes) at end of SAC file: %s (npts=%i)'
322 % (len(filedata) - nbh+nblocks*nbb, filename, self.npts))
324 def write(self, filename, byte_sex='little'):
325 '''
326 Write SAC file.
327 '''
329 self.check()
331 # create header data
332 format = SacFile.header_num_format[byte_sex]
333 numerical_values = []
334 string_values = []
335 for k in SacFile.header_keys:
336 v = self.__dict__[k]
337 vv = self.val_for_file(k, v)
338 if SacFile.t_lookup[k] == 'k':
339 string_values.append(vv)
340 else:
341 numerical_values.append(vv)
343 header_data = struct.pack(format, *numerical_values)
344 header_data += bytes(''.join(string_values).encode('ascii'))
346 # check that enough data is available
347 nblocks = self.ndatablocks()
348 if len(self.data) != nblocks:
349 raise SacError(
350 'Need %i data blocks for file type %s.'
351 % (nblocks, SacFile.header_num2name[self.iftype]))
353 for fdata in self.data:
354 if len(fdata) != self.npts:
355 raise SacError(
356 'Data length (%i) does not match NPTS header value (%i)'
357 % (len(fdata), self.npts))
359 # dump data to file
360 with open(filename, 'wb') as f:
361 f.write(header_data)
362 for fdata in self.data:
363 f.write(fdata.astype(num.float32).tostring())
365 def __str__(self):
366 str = ''
367 for k in SacFile.header_keys:
368 v = self.__dict__[k]
369 if v is not None:
370 if k in SacFile.enum_header_vars:
371 if v in SacFile.header_num2name:
372 v = SacFile.header_num2name[v]
373 str += '%s: %s\n' % (k, v)
375 return str
377 def to_trace(self):
379 assert self.iftype == SacFile.header_name2num['itime']
380 assert self.leven
382 tmin = self.get_ref_time() + self.b
383 tmax = tmin + self.delta*(self.npts-1)
385 data = None
386 if self.data:
387 data = self.data[0]
389 meta = {}
390 exclude = ('b', 'e', 'knetwk', 'kstnm', 'khole', 'kcmpnm', 'delta',
391 'nzyear', 'nzjday', 'nzhour', 'nzmin', 'nzsec', 'nzmsec')
393 for k in SacFile.header_keys:
394 if k in exclude:
395 continue
396 v = self.__dict__[k]
397 if v is not None:
398 meta[reuse(k)] = v
400 return trace.Trace(
401 nonetoempty(self.knetwk)[:2],
402 nonetoempty(self.kstnm)[:5],
403 nonetoempty(self.khole)[:2],
404 nonetoempty(self.kcmpnm)[:3],
405 tmin,
406 tmax,
407 self.delta,
408 data,
409 meta=meta)
412def iload(filename, load_data=True):
414 try:
415 sacf = SacFile(filename, load_data=load_data)
416 tr = sacf.to_trace()
417 yield tr
419 except (OSError, SacError) as e:
420 raise FileLoadError(e)
423def detect(first512):
425 if len(first512) < 512: # SAC header is 632 bytes long
426 return False
428 for sex in 'little', 'big':
429 format = SacFile.header_num_format[sex]
430 nbn = struct.calcsize(format)
432 hv = list(struct.unpack(format, first512[:nbn]))
433 iftype, nvhdr, npts, leven, delta, e, b = [
434 hv[i] for i in (85, 76, 79, 105, 0, 6, 5)]
436 if iftype not in [SacFile.header_name2num[x] for x in (
437 'itime', 'irlim', 'iamph', 'ixy', 'ixyz')]:
438 continue
439 if nvhdr < 1 or 20 < nvhdr:
440 continue
441 if npts < 0:
442 continue
443 if leven not in (0, 1, -12345):
444 continue
445 if leven and delta <= 0.0:
446 continue
447 if e != -12345.0 and b > e:
448 continue
450 return True
452 return False