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