Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/segy.py: 80%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Reader for SEGY files. 

8''' 

9 

10import numpy as num 

11import struct 

12import calendar 

13 

14from pyrocko import trace, util 

15from .io_common import FileLoadError 

16 

17 

18def ibm2ieee(ibm): 

19 """ 

20 Converts an IBM floating point number into IEEE format. 

21 :param: ibm - 32 bit unsigned integer: unpack('>L', f.read(4)) 

22 """ 

23 if ibm == 0: 

24 return 0.0 

25 sign = ibm >> 31 & 0x01 

26 exponent = ibm >> 24 & 0x7f 

27 mantissa = (ibm & 0x00ffffff) / float(pow(2, 24)) 

28 print('x', sign, exponent - 64, mantissa) 

29 return (1 - 2 * sign) * mantissa * float(pow(16, exponent - 64)) 

30 

31 

32def unpack_ibm_f4(data): 

33 ibm = num.frombuffer(data, dtype='>u4').astype(num.int32) 

34 sign = (ibm >> 31) & 0x01 

35 exponent = (ibm >> 24) & 0x7f 

36 mantissa = (ibm & 0x00ffffff) / float(pow(2, 24)) 

37 xxx = (1 - 2 * sign) * mantissa * (16.0 ** (exponent - 64))\ 

38 .astype(float) 

39 # for i in range(len(data)/4): 

40 # yyy = ibm2ieee(struct.unpack('>L', data[i*4:(i+1)*4])[0]) 

41 # print('y', sign[i], exponent[i] - 64, mantissa[i]) 

42 # print(xxx[i], yyy) 

43 # if xxx[i] != yyy: 

44 # sys.exit() 

45 # print 

46 return xxx 

47 

48 

49class SEGYError(Exception): 

50 pass 

51 

52 

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

54 ''' 

55 Read SEGY file. 

56 

57 filename -- Name of SEGY file. 

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

59 ''' 

60 

61 endianness = endianness 

62 

63 nbth = 3200 

64 nbbh = 400 

65 nbthx = 3200 

66 nbtrh = 240 

67 

68 try: 

69 f = open(filename, 'rb') 

70 

71 textual_file_header = f.read(nbth) 

72 

73 if len(textual_file_header) != nbth: 

74 raise SEGYError('incomplete textual file header') 

75 

76 binary_file_header = f.read(nbbh) 

77 if len(binary_file_header) != nbbh: 

78 raise SEGYError('incomplete binary file header') 

79 

80 line_number = struct.unpack(endianness+'1I', binary_file_header[4:8]) 

81 hvals = struct.unpack(endianness+'24H', binary_file_header[12:12+24*2]) 

82 (ntraces, nauxtraces, deltat_us, deltat_us_orig, nsamples, 

83 nsamples_orig, format, ensemble_fold) = hvals[0:8] 

84 

85 (segy_revision, fixed_length_traces, nextended_headers) = \ 

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

87 

88 if ntraces == 0 and nauxtraces == 0: 

89 ntraces = 1 

90 

91 formats = { 

92 1: (unpack_ibm_f4, 4, '4-byte IBM floating-point'), 

93 2: (endianness+'i4', 4, "4-byte, two's complement integer"), 

94 3: (endianness+'i4', 2, "2-byte, two's complement integer"), 

95 4: (None, 4, '4-byte fixed-point with gain (obolete)'), 

96 5: (endianness+'f4', 4, '4-byte IEEE floating-point'), 

97 6: (None, 0, 'not currently used'), 

98 7: (None, 0, 'not currently used'), 

99 8: ('i1', 1, "1-byte, two's complement integer")} 

100 

101 dtype = formats[format][0] 

102 sample_size = formats[format][1] 

103 if dtype is None: 

104 raise SEGYError('unsupported sample data format %i: %s' % ( 

105 format, formats[format][2])) 

106 

107 for ihead in range(nextended_headers): 

108 f.read(nbthx) 

109 

110 atend = False 

111 while not atend: 

112 for itrace in range((ntraces+nauxtraces)): 

113 trace_header = f.read(nbtrh) 

114 if len(trace_header) == 0: 

115 atend = True 

116 break 

117 

118 if len(trace_header) != nbtrh: 

119 raise SEGYError('incomplete trace header') 

120 

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

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

123 

124 (ensemblex, ensembley) = \ 

125 struct.unpack(endianness+'2f', trace_header[180:180+2*4]) 

126 (ensemble_num,) = \ 

127 struct.unpack(endianness+'1I', trace_header[20:24]) 

128 (trensemble_num,) = \ 

129 struct.unpack(endianness+'1I', trace_header[24:28]) 

130 (trace_number,) = \ 

131 struct.unpack(endianness+'1I', trace_header[0:4]) 

132 (trace_numbersegy,) = \ 

133 struct.unpack(endianness+'1I', trace_header[4:8]) 

134 (orfield_num,) = \ 

135 struct.unpack(endianness+'1I', trace_header[8:12]) 

136 (ortrace_num,) = \ 

137 struct.unpack(endianness+'1I', trace_header[12:16]) 

138 

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

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

141 

142 tvals = struct.unpack( 

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

144 

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

146 

147 tscalar = struct.unpack( 

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

149 

150 if tscalar == 0: 

151 tscalar = 1. 

152 elif tscalar < 0: 

153 tscalar = 1.0 / tscalar 

154 else: 

155 tscalar = float(tscalar) 

156 

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

158 

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

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

161 

162 # maybe not standard? 

163 (msecs, usecs) = struct.unpack( 

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

165 

166 try: 

167 if year < 100: 

168 year += 2000 

169 

170 tmin = util.to_time_float(calendar.timegm( 

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

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

173 

174 except Exception: 

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

176 

177 if fixed_length_traces: 

178 if (nsamples_this, deltat_us_this) \ 

179 != (nsamples, deltat_us): 

180 

181 raise SEGYError( 

182 'trace of incorrect length or sampling ' 

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

184 

185 if load_data: 

186 datablock = f.read(nsamples_this*sample_size) 

187 if len(datablock) != nsamples_this*sample_size: 

188 raise SEGYError('incomplete trace data') 

189 

190 if isinstance(dtype, str): 

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

192 else: 

193 data = dtype(datablock) 

194 

195 tmax = None 

196 else: 

197 f.seek(nsamples_this*sample_size, 1) 

198 tmax = tmin + deltat_us_this/1000000.*(nsamples_this-1) 

199 data = None 

200 

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

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

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

204 

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

206 

207 tr = trace.Trace( 

208 '', 

209 '%05i' % (line_number), 

210 '%02i' % (ensemble_num), 

211 '%03i' % (ortrace_num), 

212 tmin=tmin, 

213 tmax=tmax, 

214 deltat=deltat_us_this/1000000., 

215 ydata=data, 

216 meta=dict( 

217 orfield_num=orfield_num, 

218 distance=float(idist))) 

219 

220 yield tr 

221 

222 except (OSError, SEGYError) as e: 

223 raise FileLoadError(e) 

224 

225 finally: 

226 f.close()