1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import division, absolute_import 

6 

7import numpy as num 

8import struct 

9import calendar 

10 

11from pyrocko import trace, util 

12from .io_common import FileLoadError 

13 

14 

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)) 

27 

28 

29def unpack_ibm_f4(data): 

30 ibm = num.frombuffer(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 

44 

45 

46class SEGYError(Exception): 

47 pass 

48 

49 

50def iload(filename, load_data, endianness='>'): 

51 ''' 

52 Read SEGY file. 

53 

54 filename -- Name of SEGY file. 

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

56 ''' 

57 

58 endianness = endianness 

59 

60 nbth = 3200 

61 nbbh = 400 

62 nbthx = 3200 

63 nbtrh = 240 

64 

65 try: 

66 f = open(filename, 'rb') 

67 

68 textual_file_header = f.read(nbth) 

69 

70 if len(textual_file_header) != nbth: 

71 raise SEGYError('incomplete textual file header') 

72 

73 binary_file_header = f.read(nbbh) 

74 if len(binary_file_header) != nbbh: 

75 raise SEGYError('incomplete binary file header') 

76 

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] 

81 

82 (segy_revision, fixed_length_traces, nextended_headers) = \ 

83 struct.unpack(endianness+'3H', binary_file_header[100:100+3*2]) 

84 

85 if ntraces == 0 and nauxtraces == 0: 

86 ntraces = 1 

87 

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")} 

97 

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])) 

103 

104 for ihead in range(nextended_headers): 

105 f.read(nbthx) 

106 

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 

114 

115 if len(trace_header) != nbtrh: 

116 raise SEGYError('incomplete trace header') 

117 

118 (scoordx, scoordy, gcoordx, gcoordy) = \ 

119 struct.unpack(endianness+'4f', trace_header[72:72+4*4]) 

120 

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]) 

135 

136 # don't know if this is standard: distance to shot [m] as int 

137 (idist,) = struct.unpack(endianness+'1I', trace_header[36:40]) 

138 

139 tvals = struct.unpack( 

140 endianness+'12H', trace_header[94:94+12*2]) 

141 

142 (nsamples_this, deltat_us_this) = tvals[-2:] 

143 

144 tscalar = struct.unpack( 

145 endianness+'1H', trace_header[214:216])[0] 

146 

147 if tscalar == 0: 

148 tscalar = 1. 

149 elif tscalar < 0: 

150 tscalar = 1.0 / tscalar 

151 else: 

152 tscalar = float(tscalar) 

153 

154 tvals = [x * tscalar for x in tvals[:-2]] 

155 

156 (year, doy, hour, minute, second) = \ 

157 struct.unpack(endianness+'5H', trace_header[156:156+2*5]) 

158 

159 # maybe not standard? 

160 (msecs, usecs) = struct.unpack( 

161 endianness+'2H', trace_header[168:168+4]) 

162 

163 try: 

164 if year < 100: 

165 year += 2000 

166 

167 tmin = util.to_time_float(calendar.timegm( 

168 (year, 1, doy, hour, minute, second))) \ 

169 + msecs * 1.0e-3 + usecs * 1.0e-6 

170 

171 except Exception: 

172 raise SEGYError('invalid start date/time') 

173 

174 if fixed_length_traces: 

175 if (nsamples_this, deltat_us_this) \ 

176 != (nsamples, deltat_us): 

177 

178 raise SEGYError( 

179 'trace of incorrect length or sampling ' 

180 'rate (trace=%i)' % itrace+1) 

181 

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') 

186 

187 if isinstance(dtype, str): 

188 data = num.frombuffer(datablock, dtype=dtype) 

189 else: 

190 data = dtype(datablock) 

191 

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 

197 

198 if data is not None and isinstance(dtype, str) and ( 

199 str(dtype).startswith('<') 

200 or str(dtype).startswith('>')): 

201 

202 data = data.astype(str(dtype)[1:]) 

203 

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))) 

216 

217 yield tr 

218 

219 except (OSError, SEGYError) as e: 

220 raise FileLoadError(e) 

221 

222 finally: 

223 f.close()