1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
7import struct
8import calendar
10from pyrocko import trace, util
11from .io_common import FileLoadError
14def ibm2ieee(ibm):
15 """
16 Converts an IBM floating point number into IEEE format.
17 :param: ibm - 32 bit unsigned integer: unpack('>L', f.read(4))
18 """
19 if ibm == 0:
20 return 0.0
21 sign = ibm >> 31 & 0x01
22 exponent = ibm >> 24 & 0x7f
23 mantissa = (ibm & 0x00ffffff) / float(pow(2, 24))
24 print('x', sign, exponent - 64, mantissa)
25 return (1 - 2 * sign) * mantissa * float(pow(16, exponent - 64))
28def unpack_ibm_f4(data):
29 ibm = num.frombuffer(data, dtype='>u4').astype(num.int32)
30 sign = (ibm >> 31) & 0x01
31 exponent = (ibm >> 24) & 0x7f
32 mantissa = (ibm & 0x00ffffff) / float(pow(2, 24))
33 xxx = (1 - 2 * sign) * mantissa * (16.0 ** (exponent - 64))\
34 .astype(float)
35 # for i in range(len(data)/4):
36 # yyy = ibm2ieee(struct.unpack('>L', data[i*4:(i+1)*4])[0])
37 # print('y', sign[i], exponent[i] - 64, mantissa[i])
38 # print(xxx[i], yyy)
39 # if xxx[i] != yyy:
40 # sys.exit()
41 # print
42 return xxx
45class SEGYError(Exception):
46 pass
49def iload(filename, load_data, endianness='>'):
50 '''
51 Read SEGY file.
53 filename -- Name of SEGY file.
54 load_data -- If True, the data is read, otherwise only read headers.
55 '''
57 endianness = endianness
59 nbth = 3200
60 nbbh = 400
61 nbthx = 3200
62 nbtrh = 240
64 try:
65 f = open(filename, 'rb')
67 textual_file_header = f.read(nbth)
69 if len(textual_file_header) != nbth:
70 raise SEGYError('incomplete textual file header')
72 binary_file_header = f.read(nbbh)
73 if len(binary_file_header) != nbbh:
74 raise SEGYError('incomplete binary file header')
76 line_number = struct.unpack(endianness+'1I', binary_file_header[4:8])
77 hvals = struct.unpack(endianness+'24H', binary_file_header[12:12+24*2])
78 (ntraces, nauxtraces, deltat_us, deltat_us_orig, nsamples,
79 nsamples_orig, format, ensemble_fold) = hvals[0:8]
81 (segy_revision, fixed_length_traces, nextended_headers) = \
82 struct.unpack(endianness+'3H', binary_file_header[100:100+3*2])
84 if ntraces == 0 and nauxtraces == 0:
85 ntraces = 1
87 formats = {
88 1: (unpack_ibm_f4, 4, "4-byte IBM floating-point"),
89 2: (endianness+'i4', 4, "4-byte, two's complement integer"),
90 3: (endianness+'i4', 2, "2-byte, two's complement integer"),
91 4: (None, 4, "4-byte fixed-point with gain (obolete)"),
92 5: (endianness+'f4', 4, "4-byte IEEE floating-point"),
93 6: (None, 0, "not currently used"),
94 7: (None, 0, "not currently used"),
95 8: ('i1', 1, "1-byte, two's complement integer")}
97 dtype = formats[format][0]
98 sample_size = formats[format][1]
99 if dtype is None:
100 raise SEGYError('unsupported sample data format %i: %s' % (
101 format, formats[format][2]))
103 for ihead in range(nextended_headers):
104 f.read(nbthx)
106 atend = False
107 while not atend:
108 for itrace in range((ntraces+nauxtraces)):
109 trace_header = f.read(nbtrh)
110 if len(trace_header) == 0:
111 atend = True
112 break
114 if len(trace_header) != nbtrh:
115 raise SEGYError('incomplete trace header')
117 (scoordx, scoordy, gcoordx, gcoordy) = \
118 struct.unpack(endianness+'4f', trace_header[72:72+4*4])
120 (ensemblex, ensembley) = \
121 struct.unpack(endianness+'2f', trace_header[180:180+2*4])
122 (ensemble_num,) = \
123 struct.unpack(endianness+'1I', trace_header[20:24])
124 (trensemble_num,) = \
125 struct.unpack(endianness+'1I', trace_header[24:28])
126 (trace_number,) = \
127 struct.unpack(endianness+'1I', trace_header[0:4])
128 (trace_numbersegy,) = \
129 struct.unpack(endianness+'1I', trace_header[4:8])
130 (orfield_num,) = \
131 struct.unpack(endianness+'1I', trace_header[8:12])
132 (ortrace_num,) = \
133 struct.unpack(endianness+'1I', trace_header[12:16])
135 # don't know if this is standard: distance to shot [m] as int
136 (idist,) = struct.unpack(endianness+'1I', trace_header[36:40])
138 tvals = struct.unpack(
139 endianness+'12H', trace_header[94:94+12*2])
141 (nsamples_this, deltat_us_this) = tvals[-2:]
143 tscalar = struct.unpack(
144 endianness+'1H', trace_header[214:216])[0]
146 if tscalar == 0:
147 tscalar = 1.
148 elif tscalar < 0:
149 tscalar = 1.0 / tscalar
150 else:
151 tscalar = float(tscalar)
153 tvals = [x * tscalar for x in tvals[:-2]]
155 (year, doy, hour, minute, second) = \
156 struct.unpack(endianness+'5H', trace_header[156:156+2*5])
158 # maybe not standard?
159 (msecs, usecs) = struct.unpack(
160 endianness+'2H', trace_header[168:168+4])
162 try:
163 if year < 100:
164 year += 2000
166 tmin = util.to_time_float(calendar.timegm(
167 (year, 1, doy, hour, minute, second))) \
168 + msecs * 1.0e-3 + usecs * 1.0e-6
170 except Exception:
171 raise SEGYError('invalid start date/time')
173 if fixed_length_traces:
174 if (nsamples_this, deltat_us_this) \
175 != (nsamples, deltat_us):
177 raise SEGYError(
178 'trace of incorrect length or sampling '
179 'rate (trace=%i)' % itrace+1)
181 if load_data:
182 datablock = f.read(nsamples_this*sample_size)
183 if len(datablock) != nsamples_this*sample_size:
184 raise SEGYError('incomplete trace data')
186 if isinstance(dtype, str):
187 data = num.frombuffer(datablock, dtype=dtype)
188 else:
189 data = dtype(datablock)
191 tmax = None
192 else:
193 f.seek(nsamples_this*sample_size, 1)
194 tmax = tmin + deltat_us_this/1000000.*(nsamples_this-1)
195 data = None
197 if data is not None and isinstance(dtype, str) and (
198 str(dtype).startswith('<')
199 or str(dtype).startswith('>')):
201 data = data.astype(str(dtype)[1:])
203 tr = trace.Trace(
204 '',
205 '%05i' % (line_number),
206 '%02i' % (ensemble_num),
207 '%03i' % (ortrace_num),
208 tmin=tmin,
209 tmax=tmax,
210 deltat=deltat_us_this/1000000.,
211 ydata=data,
212 meta=dict(
213 orfield_num=orfield_num,
214 distance=float(idist)))
216 yield tr
218 except (OSError, SEGYError) as e:
219 raise FileLoadError(e)
221 finally:
222 f.close()