SAC Tutorial and Case Study

headersac.csv

Word	Type	NAMES	o	o	o	o
0	F	DELTA	DEPMIN	DEPMAX	SCALE	ODELTA
5	F	B	E	O	A	INTERNAL
10	F	T0	T1	T2	T3	T4
15	F	T5	T6	T7	T8	T9
20	F	F	RESP0	RESP1	RESP2	RESP3
25	F	RESP4	RESP5	RESP6	RESP7	RESP8
30	F	RESP9	STLA	STLO	STEL	STDP
35	F	EVLA	EVLO	EVEL	EVDP	MAG
40	F	USER0	USER1	USER2	USER3	USER4
45	F	USER5	USER6	USER7	USER8	USER9
50	F	DIST	AZ	BAZ	GCARC	INTERNAL
55	F	INTERNAL	DEPMEN	CMPAZ	CMPINC	XMINIMUM
60	F	XMAXIMUM	YMINIMUM	YMAXIMUM	UNUSED	UNUSED
65	F	UNUSED	UNUSED	UNUSED	UNUSED	UNUSED
70	I	NZYEAR	NZJDAY	NZHOUR	NZMIN	NZSEC
75	I	NZMSEC	NVHDR	NORID	NEVID	NPTS
80	I	INTERNAL	NWFID	NXSIZE	NYSIZE	UNUSED
85	I	IFTYPE	IDEP	IZTYPE	UNUSED	IINST
90	I	ISTREG	IEVREG	IEVTYP	IQUAL	ISYNTH
95	I	IMAGTYP	IMAGSRC	UNUSED	UNUSED	UNUSED
100	I	UNUSED	UNUSED	UNUSED	UNUSED	UNUSED
105	L	LEVEN	LPSPOL	LOVROK	LCALDA	UNUSED
110	K	KSTNM	KEVNM*	NaN	NaN	NaN
116	K	KHOLE	KO	KA	NaN	NaN
122	K	KT0	KT1	KT2	NaN	NaN
128	K	KT3	KT4	KT5	NaN	NaN
134	K	KT6	KT7	KT8	NaN	NaN
140	K	KT9	KF	KUSER0	NaN	NaN
146	K	KUSER1	KUSER2	KCMPNM	NaN	NaN
152	K	KNETWK	KDATRD	KINST	NaN	NaN

readbinarysac.py

# import os, struct, numpy and matplotlib
import os
import struct
import numpy as np
import matplotlib.pyplot as plt

# identify header sac filename and a sac binary filename
headerfile = "headersac.csv"
filename = "example.sac"

# read header sac filename and change to the list python
file = open(headerfile,'r')
hdr = file.readlines()
for i in range(len(hdr)):
    hdr[i] = hdr[i].split()
file.close()

# clustering and making the categories the header type
# hdrVar as a list of SAC header type items
# spec as a list of struct type corresponded and length of characters
hdrVar = []
spec = []
for i in range(len(hdr)):
    for j in range(len(hdr[0])):
        if i != 0:
            if j != 0 and j != 1:
                if hdr[i][j] != 'NaN':
                    hdrVar.append(hdr[i][j])
                    if hdr[i][1] == 'F':
                        spec.append(['f',4])
                    if hdr[i][1] == 'I':
                        spec.append(['i',4])
                    if hdr[i][1] == 'L':
                        spec.append(['?',4])
                    if hdr[i][1] == 'K' and hdr[i][j] != 'KEVNM*':
                        spec.append(['s',8])
                    if hdr[i][1] == 'K' and hdr[i][j] == 'KEVNM*':
                        spec.append(['s', 16])

# opening the sac binary file and storing the length to variable 'size'
f = open(filename, 'rb')
f.seek(0, os.SEEK_END)
size = f.tell()
print('size of the data is ',size)

# opening again and converting binary to format characters
# bi variable as a control of what the length that has been read
# amp variable as a list of amplitude data
# the list of SAC headers only display in the command/output window
# except DELTA due to plotting
f = open(filename,'rb')
bi = 0
amp = []
DELTA = 0
while bi < size:
    if bi < len(hdrVar):
        for i in range(len(hdrVar)):
            if spec[i][0] != 's' and spec[i][0] != '?':
                datBin = f.read(spec[i][1])
                print(hdrVar[i], struct.unpack('<{}'.format(spec[i][0]),datBin)[0])
                if hdrVar[i] == 'DELTA':
                    DELTA = struct.unpack('<{}'.format(spec[i][0]),datBin)[0]
                bi += spec[i][1]
            elif spec[i][0] == '?':
                datBin = f.read(spec[i][1])
                print(hdrVar[i], struct.unpack('<{}'.format(spec[i][0]), datBin[:1])[0])
                bi += spec[i][1]
            else:
                datBin = f.read(spec[i][1])
                print(hdrVar[i], str(struct.unpack('<{}{}'.format(spec[i][1],spec[i][0]), datBin)[0]).split("'")[1])
                bi += spec[i][1]
    else:
        datBin = f.read(4)
        amp.append(struct.unpack('<f', datBin))
        bi += 4

# plot the data
# t as the time
t = np.r_[0:len(amp)*DELTA:DELTA]
plt.plot(t,amp)
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.title('Example of plot data')
plt.ylabel('Amplitude')
plt.xlabel('Time (s)')
plt.show()