Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/elastic10_to_rotational8.py: 0%

1import os 

2import logging 

3import shutil 

4from pyrocko import gf, guts, util 

5import numpy as num 

6 

7 

8logger = logging.getLogger('main') 

9 

10 

11def a(*args): 

12 return num.array(args, dtype=float) 

13 

14 

15def ai(*args): 

16 return num.array(args, dtype=int) 

17 

18 

19def almost_equal(a, b, eps): 

20 return abs(a - b < eps) 

21 

22 

23def configure_diff(x, xmin, xmax, dx): 

24 if almost_equal(x, xmin, 1e-6*dx): 

25 return dx, 0. 

26 elif almost_equal(x, xmax, 1e-6*dx): 

27 return 0., dx 

28 else: 

29 return dx, dx 

30 

31 

32def sindex(args): 

33 return '(%s)' % ', '.join('%g' % x for x in args) 

34 

35 

36def elastic10_to_rotational8( 

37 source_store_dir, 

38 dest_store_dir, 

39 show_progress=True): 

40 

41 def _raise(message): 

42 raise gf.StoreError('elastic10_to_rotational8: %s' % message) 

43 

44 source = gf.Store(source_store_dir) 

45 

46 if source.config.component_scheme not in ('elastic10', 'elastic10_fd'): 

47 _raise( 

48 'Only `elastic10` and `elastic10_fd` component schemes supported ' 

49 'for input store.') 

50 

51 have_fd_components = source.config.component_scheme == 'elastic10_fd' 

52 

53 if source.config.effective_stored_quantity != 'displacement': 

54 _raise('Stored quantity of input store must be displacement.') 

55 

56 if not os.path.exists(dest_store_dir): 

57 config = guts.clone(source.config) 

58 config.id += '_rotfd' 

59 config.__dict__['ncomponents'] = 8 

60 config.__dict__['component_scheme'] = 'rotational8' 

61 config.stored_quantity = 'rotation_displacement' 

62 config.fd_distance_delta = None 

63 if config.short_type == 'B': 

64 config.fd_receiver_depth_delta = None 

65 

66 util.ensuredirs(dest_store_dir) 

67 gf.Store.create(dest_store_dir, config=config) 

68 shutil.copytree( 

69 os.path.join(source_store_dir, 'phases'), 

70 os.path.join(dest_store_dir, 'phases')) 

71 

72 try: 

73 gf.store.Store.create_dependants(dest_store_dir) 

74 except gf.StoreError: 

75 pass 

76 

77 dest = gf.Store(dest_store_dir, 'w') 

78 

79 ta = source.config.short_type 

80 tb = dest.config.short_type 

81 

82 if have_fd_components and not (ta, tb) == ('A', 'A'): 

83 _raise( 

84 'For input stores with component scheme `elastic10_fd`, only ' 

85 'conversion from store type A to A is currently supported.') 

86 

87 if (ta, tb) not in (('A', 'A'), ('B', 'B'), ('B', 'A')): 

88 _raise('Cannot convert type %s to type %s.' % (ta, tb)) 

89 

90 if source.config.distance_delta > dest.config.distance_delta: 

91 _raise( 

92 'Distance spacing of output store must be equal to or larger ' 

93 'distance spacing of input store.') 

94 

95 if (ta, tb) == ('B', 'B'): 

96 if source.config.receiver_depth_delta \ 

97 > dest.config.receiver_depth_delta: 

98 _raise( 

99 'Depth spacing of output store must be equal to or larger ' 

100 'depth spacing of input store.') 

101 

102 if ta == 'A': 

103 if source.config.receiver_depth \ 

104 != source.config.earthmodel_1d.discontinuity('surface').z: 

105 _raise( 

106 'When input store is of type A, receivers must be at the ' 

107 'surface.') 

108 

109 if have_fd_components: 

110 dx = source.config.fd_distance_delta 

111 else: 

112 dx = source.config.distance_delta 

113 

114 if ta == 'B': 

115 dz = source.config.receiver_depth_delta 

116 else: 

117 dz = None 

118 

119 if show_progress: 

120 pbar = util.progressbar( 

121 'elastic10_to_rotational8', 

122 source.config.nrecords/source.config.ncomponents) 

123 

124 try: 

125 for i, args in enumerate(dest.config.iter_nodes(level=-1)): 

126 if ta == 'A': 

127 sz, x = args 

128 

129 if have_fd_components: 

130 dx1 = dx2 = dx 

131 else: 

132 dx1, dx2 = configure_diff( 

133 x, 

134 source.config.distance_min, 

135 source.config.distance_max, 

136 dx) 

137 

138 odx = dx1+dx2 

139 

140 for ig, (sum_args, weights) in enumerate(zip( 

141 [ 

142 (a(), a(), ai()), 

143 (a(), a(), ai()), 

144 (a(sz, sz), a(x, x), ai(25, 5)), 

145 (a(sz, sz), a(x, x), ai(26, 6)), 

146 (a(sz, sz), a(x, x), ai(27, 7)), 

147 (a(sz, sz), a(x, x), ai(23, 3)), 

148 (a(sz, sz), a(x, x), ai(24, 4)), 

149 (a(sz, sz), a(x, x), ai(29, 9)), 

150 ] if have_fd_components else [ 

151 (a(), a(), ai()), 

152 (a(), a(), ai()), 

153 (a(sz, sz), a(x+dx1, x-dx2), ai(5, 5)), 

154 (a(sz, sz), a(x+dx1, x-dx2), ai(6, 6)), 

155 (a(sz, sz), a(x+dx1, x-dx2), ai(7, 7)), 

156 (a(sz, sz), a(x+dx1, x-dx2), ai(3, 3)), 

157 (a(sz, sz), a(x+dx1, x-dx2), ai(4, 4)), 

158 (a(sz, sz), a(x+dx1, x-dx2), ai(9, 9)), 

159 ], 

160 [ 

161 a(), 

162 a(), 

163 1.0 / a(-odx, odx), 

164 1.0 / a(-odx, odx), 

165 1.0 / a(-odx, odx), 

166 0.5 / a(odx, -odx), 

167 0.5 / a(odx, -odx), 

168 1.0 / a(-odx, odx), 

169 ])): 

170 

171 if weights.size == 0: 

172 h = gf.GFTrace(is_zero=True, itmin=0) 

173 else: 

174 h = source.sum( 

175 sum_args, 

176 num.zeros_like(weights), 

177 weights, 

178 optimization='disable') 

179 

180 dest.put(args + (ig,), h) 

181 else: 

182 rz, sz, x = args 

183 

184 dx1, dx2 = configure_diff( 

185 x, 

186 source.config.distance_min, 

187 source.config.distance_max, 

188 dx) 

189 

190 dz1, dz2 = configure_diff( 

191 rz, 

192 source.config.receiver_depth_min, 

193 source.config.receiver_depth_max, 

194 dz) 

195 

196 odx = dx1+dx2 

197 odz = dz1+dz2 

198 

199 for ig, (sum_args, weights) in enumerate(zip( 

200 [ 

201 (a(rz+dz1, rz-dz2), a(sz, sz), a(x, x), ai(3, 3)), 

202 (a(rz+dz1, rz-dz2), a(sz, sz), a(x, x), ai(4, 4)), 

203 (a(rz+dz1, rz-dz2, rz, rz), a(sz, sz, sz, sz), a(x, x, x+dx1, x-dx2), ai(0, 0, 5, 5)), # noqa 

204 (a(rz+dz1, rz-dz2, rz, rz), a(sz, sz, sz, sz), a(x, x, x+dx1, x-dx2), ai(1, 1, 6, 6)), # noqa 

205 (a(rz+dz1, rz-dz2, rz, rz), a(sz, sz, sz, sz), a(x, x, x+dx1, x-dx2), ai(2, 2, 7, 7)), # noqa 

206 (a(rz, rz), a(sz, sz), a(x+dx1, x-dx2), ai(3, 3)), 

207 (a(rz, rz), a(sz, sz), a(x+dx1, x-dx2), ai(4, 4)), 

208 (a(rz+dz1, rz-dz2, rz, rz), a(sz, sz, sz, sz), a(x, x, x+dx1, x-dx2), ai(8, 8, 9, 9)), # noqa 

209 ], 

210 [ 

211 0.5 / a(-odz, odz), 

212 0.5 / a(-odz, odz), 

213 0.5 / a(odz, -odz, -odx, odx), 

214 0.5 / a(odz, -odz, -odx, odx), 

215 0.5 / a(odz, -odz, -odx, odx), 

216 0.5 / a(odx, -odx), 

217 0.5 / a(odx, -odx), 

218 0.5 / a(odz, -odz, -odx, odx), 

219 ])): 

220 

221 h = source.sum( 

222 sum_args, 

223 num.zeros_like(weights), 

224 weights, 

225 optimization='disable') 

226 

227 dest.put(args + (ig,), h) 

228 

229 if show_progress: 

230 pbar.update(i+1) 

231 

232 finally: 

233 if show_progress: 

234 pbar.finish()