1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

# http://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

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

'''SAC IO library for Python''' 

from __future__ import absolute_import 

from builtins import range 

 

import struct 

import logging 

import math 

import numpy as num 

 

from calendar import timegm 

from time import gmtime 

 

from pyrocko import trace 

from pyrocko.util import reuse 

from .io_common import FileLoadError 

 

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

 

 

def fixdoublefloat(x): 

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

return round(x/f)*f 

 

 

class SacError(Exception): 

pass 

 

 

def nonetoempty(s): 

if s is None: 

return '' 

else: 

return s.strip() 

 

 

class SacFile(object): 

nbytes_header = 632 

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

 

header_keys = ''' 

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

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

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

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

xmaximum yminimum ymaximum unused0 unused1 unused2 unused3 unused4 unused5 

unused6 nzyear nzjday nzhour nzmin nzsec nzmsec nvhdr norid nevid npts 

internal3 nwfid nxsize nysize unused7 iftype idep iztype unused8 iinst istreg 

ievreg ievtyp iqual isynth imagtyp imagsrc unused9 unused10 unused11 unused12 

unused13 unused14 unused15 unused16 leven lpspol lovrok lcalda unused17 kstnm 

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

kuser2 kcmpnm knetwk kdatrd kinst 

'''.split() 

 

header_enum_symbols = ''' 

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

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

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

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

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

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

'''.split() 

 

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

'isynth'.split() 

 

header_num2name = dict( 

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

header_name2num = dict( 

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

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

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

ldefaults = { 

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

 

t_lookup = dict(zip(header_keys, header_types)) 

 

u_lookup = dict() 

for k in header_keys: 

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

 

def ndatablocks(self): 

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

nblocks = { 

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

}[SacFile.header_num2name[self.iftype]] 

 

if nblocks == 1 and not self.leven: 

nblocks = 2 # not sure about this... 

 

return nblocks 

 

def val_or_none(self, k, v): 

'''Replace SAC undef flags with None.''' 

if SacFile.u_lookup[k] == v: 

return None 

else: 

return v 

 

def get_ref_time(self): 

'''Get reference time as standard Unix timestamp.''' 

 

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

self.nzsec, self.nzmsec): 

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

 

return timegm((self.nzyear, 1, self.nzjday, self.nzhour, self.nzmin, 

self.nzsec)) + self.nzmsec/1000. 

 

def set_ref_time(self, timestamp): 

''' 

Set all header values defining reference time based on standard Unix 

timestamp. 

''' 

 

secs = math.floor(timestamp) 

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

if msec == 1000: 

secs += 1 

msec = 0 

 

t = gmtime(secs) 

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

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

self.nzmsec = msec 

 

def val_for_file(self, k, v): 

''' 

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

SAC file. 

''' 

 

t = SacFile.t_lookup[k] 

if v is None: 

if t == 'l': 

return SacFile.ldefaults[k] 

v = SacFile.u_lookup[k] 

if t == 'f': 

return float(v) 

elif t == 'i': 

return int(v) 

elif t == 'l': 

if v: 

return 1 

return 0 

elif t == 'k': 

ln = 8 

if k == 'kevnm': 

ln = 16 # only this header val has different length 

return v.ljust(ln)[:ln] 

 

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

if 'from_trace' in kwargs: 

self.clear() 

trace = kwargs['from_trace'] 

if trace.meta: 

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

if k in SacFile.header_keys: 

setattr(self, k, v) 

 

self.knetwk = trace.network 

self.kstnm = trace.station 

self.khole = trace.location 

self.kcmpnm = trace.channel 

self.set_ref_time(trace.tmin) 

self.delta = trace.deltat 

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

self.npts = trace.ydata.size 

self.b = 0.0 

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

 

elif args: 

self.read(*args, **kwargs) 

else: 

self.clear() 

 

def clear(self): 

'''Empty file record.''' 

 

for k in SacFile.header_keys: 

self.__dict__[k] = None 

 

# set the required attributes 

self.nvhdr = 6 

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

self.leven = True 

self.delta = 1.0 

self.npts = 0 

self.b = 0.0 

self.e = 0.0 

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

 

def check(self): 

'''Check the required header variables to have reasonable values.''' 

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

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

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

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

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

if self.npts < 0: 

raise SacError( 

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

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

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

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

raise SacError( 

'Header value DELTA should be positive for evenly spaced ' 

'samples') 

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

raise SacError( 

'Beginning value of independent variable greater than its ' 

'ending value.') 

if self.nvhdr != 6: 

logging.warn( 

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

'This file has header version %i. ' 

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

 

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

''' 

Read SAC file. 

 

filename -- Name of SAC file. 

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

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

''' 

 

nbh = SacFile.nbytes_header 

 

# read in all data 

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

if load_data: 

filedata = f.read() 

else: 

filedata = f.read(nbh) 

 

if len(filedata) < nbh: 

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

 

# possibly try out different endiannesses 

if byte_sex == 'try': 

sexes = ('little', 'big') 

else: 

sexes = (byte_sex,) 

 

for isex, sex in enumerate(sexes): 

format = SacFile.header_num_format[sex] 

nbn = struct.calcsize(format) 

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

 

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

hv.append(strings[:8].rstrip()) 

hv.append(strings[8:24].rstrip()) 

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

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

 

self.header_vals = hv 

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

vn = self.val_or_none(k, v) 

if isinstance(vn, str): 

ipos = vn.find('\x00') 

if ipos != -1: 

vn = vn[:ipos] 

 

self.__dict__[k] = vn.rstrip() 

else: 

self.__dict__[k] = vn 

 

if self.leven == -12345: 

self.leven = True 

 

self.data = [] 

try: 

self.check() 

break 

except SacError as e: 

if isex == len(sexes)-1: 

raise e 

 

self.delta = fixdoublefloat(self.delta) 

 

if byte_sex == 'try': 

logger.debug( 

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

 

# possibly get data 

if load_data: 

nblocks = self.ndatablocks() 

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

for iblock in range(nblocks): 

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

raise SacError('File is incomplete.') 

 

if sex == 'big': 

dtype = num.dtype('>f4') 

else: 

dtype = num.dtype('<f4') 

 

self.data.append(num.array(num.fromstring( 

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

dtype=dtype), 

dtype=num.float)) 

 

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

logger.warning( 

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

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

 

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

'''Write SAC file.''' 

 

self.check() 

 

# create header data 

format = SacFile.header_num_format[byte_sex] 

numerical_values = [] 

string_values = [] 

for k in SacFile.header_keys: 

v = self.__dict__[k] 

vv = self.val_for_file(k, v) 

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

string_values.append(vv) 

else: 

numerical_values.append(vv) 

 

header_data = struct.pack(format, *numerical_values) 

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

 

# check that enough data is available 

nblocks = self.ndatablocks() 

if len(self.data) != nblocks: 

raise SacError( 

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

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

 

for fdata in self.data: 

if len(fdata) != self.npts: 

raise SacError( 

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

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

 

# dump data to file 

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

f.write(header_data) 

for fdata in self.data: 

f.write(fdata.astype(num.float32).tostring()) 

 

def __str__(self): 

str = '' 

for k in SacFile.header_keys: 

v = self.__dict__[k] 

if v is not None: 

if k in SacFile.enum_header_vars: 

if v in SacFile.header_num2name: 

v = SacFile.header_num2name[v] 

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

 

return str 

 

def to_trace(self): 

 

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

assert self.leven 

 

tmin = self.get_ref_time() + self.b 

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

 

data = None 

if self.data: 

data = self.data[0] 

 

meta = {} 

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

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

 

for k in SacFile.header_keys: 

if k in exclude: 

continue 

v = self.__dict__[k] 

if v is not None: 

meta[reuse(k)] = v 

 

return trace.Trace( 

nonetoempty(self.knetwk)[:2], 

nonetoempty(self.kstnm)[:5], 

nonetoempty(self.khole)[:2], 

nonetoempty(self.kcmpnm)[:3], 

tmin, 

tmax, 

self.delta, 

data, 

meta=meta) 

 

 

def iload(filename, load_data=True): 

 

try: 

sacf = SacFile(filename, load_data=load_data) 

tr = sacf.to_trace() 

yield tr 

 

except (OSError, SacError) as e: 

raise FileLoadError(e) 

 

 

def detect(first512): 

 

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

return False 

 

for sex in 'little', 'big': 

format = SacFile.header_num_format[sex] 

nbn = struct.calcsize(format) 

 

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

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

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

 

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

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

continue 

if nvhdr < 1 or 20 < nvhdr: 

continue 

if npts < 0: 

continue 

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

continue 

if leven and delta <= 0.0: 

continue 

if e != -12345.0 and b > e: 

continue 

 

return True 

 

return False