Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/phase_pick/plot.py: 13%

179 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-26 16:25 +0000

1import math 

2 

3from matplotlib import pyplot as plt, colors as mcolors 

4import numpy as num 

5 

6from pyrocko import util 

7from pyrocko.guts import Tuple, Float 

8from pyrocko.plot import mpl_init, mpl_margins 

9 

10from grond import meta 

11from grond.plot.config import PlotConfig 

12from grond.plot.collection import PlotItem 

13 

14from .target import PhasePickTarget 

15from ..plot import StationDistributionPlot 

16 

17km = 1000. 

18d2r = math.pi / 180. 

19 

20 

21class PickResidualsPlot(PlotConfig): 

22 '''Travel time residuals plot.''' 

23 

24 name = 'pick_residuals' 

25 

26 size_cm = Tuple.T( 

27 2, Float.T(), 

28 default=(15.0, 10.0), 

29 help='Width and height of the figure in [cm]') 

30 

31 def make(self, environ): 

32 environ.setup_modelling() 

33 

34 cm = environ.get_plot_collection_manager() 

35 mpl_init(fontsize=self.font_size) 

36 

37 ds = environ.get_dataset() 

38 history = environ.get_history(subset='harvest') 

39 optimiser = environ.get_optimiser() 

40 

41 cm.create_group_mpl( 

42 self, 

43 self.draw_figures(ds, history, optimiser), 

44 title=u'Pick Residuals', 

45 section='fits', 

46 feather_icon='watch', 

47 description=u''' 

48Travel time residuals. 

49 

50The difference between observed and synthetic phase arrival times is shown 

51for a subset of the models in the solution ensemble. Models with good global 

52performance are shown in red, less good ones in blue. Vertical black bars mark 

53the classical best model. 

54''') 

55 

56 def draw_figures(self, ds, history, optimiser): 

57 show_mean_residuals = False 

58 

59 # gms = history.get_sorted_primary_misfits()[::-1] 

60 models = history.get_sorted_primary_models()[::-1] 

61 problem = history.problem 

62 

63 targets = [ 

64 t for t in problem.targets if 

65 isinstance(t, PhasePickTarget)] 

66 

67 # targets.sort(key=lambda t: t.distance_to(problem.base_source)) 

68 # tpaths = sorted(set(t.path for t in targets)) 

69 

70 ntargets = len(targets) 

71 nmodels = history.nmodels 

72 tts = num.zeros((nmodels, ntargets, 2)) 

73 distances = num.zeros((nmodels, ntargets)) 

74 azimuths = num.zeros((nmodels, ntargets)) 

75 

76 for imodel in range(nmodels): 

77 model = models[imodel, :] 

78 source = problem.get_source(model) 

79 results = problem.evaluate(model, targets=targets) 

80 for itarget, result in enumerate(results): 

81 result = results[itarget] 

82 tts[imodel, itarget, :] = result.tobs, result.tsyn 

83 distances[imodel, itarget] = source.distance_to( 

84 targets[itarget]) 

85 

86 azimuths[imodel, itarget] = source.azibazi_to( 

87 targets[itarget])[0] 

88 

89 ok = num.all(num.isfinite(tts), axis=2) 

90 ok = num.all(ok, axis=0) 

91 

92 targets = [ 

93 target for (itarget, target) in enumerate(targets) if ok[itarget]] 

94 tts = tts[:, ok, :] 

95 distances = distances[:, ok] 

96 azimuths = azimuths[:, ok] 

97 residuals = tts[:, :, 0] - tts[:, :, 1] 

98 residual_amax = num.max(num.abs(residuals)) 

99 mean_residuals = num.mean(residuals, axis=0) 

100 mean_distances = num.mean(distances, axis=0) 

101 

102 isort = num.array( 

103 [x[2] for x in sorted( 

104 zip(targets, mean_distances, range(len(targets))), 

105 key=lambda x: (x[0].path, x[1]))], dtype=int) 

106 

107 distances = distances[:, isort] 

108 distance_min = num.min(distances) 

109 distance_max = num.max(distances) 

110 azimuths = azimuths[:, isort] 

111 targets = [targets[i] for i in isort] 

112 ntargets = len(targets) 

113 residuals = residuals[:, isort] 

114 mean_residuals = mean_residuals[isort] 

115 

116 icolor = num.arange(nmodels) 

117 

118 norm = mcolors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor)) 

119 cmap = plt.get_cmap('coolwarm') 

120 

121 napprox = max( 

122 1, int(round((self.size_points[1] / self.font_size - 7)))) 

123 

124 npages = (ntargets-1) // napprox + 1 

125 ntargets_page = (ntargets-1) // npages + 1 

126 

127 for ipage in range(npages): 

128 

129 ilo, ihi = ipage*ntargets_page, (ipage+1)*ntargets_page 

130 ihi = min(len(targets), ihi) 

131 

132 targets_page = targets[ilo:ihi] 

133 item = PlotItem( 

134 name='fig_%02i' % ipage) 

135 

136 item.attributes['targets'] = [t.string_id() for t in targets_page] 

137 

138 fig = plt.figure(figsize=self.size_inch) 

139 

140 labelpos = mpl_margins( 

141 fig, nw=2, nh=1, left=12., right=1., bottom=5., top=1., 

142 wspace=1., units=self.font_size) 

143 

144 axes1 = fig.add_subplot(1, 3, 1) 

145 axes2 = fig.add_subplot(1, 3, 2) 

146 axes3 = fig.add_subplot(1, 3, 3) 

147 

148 for axes in (axes1, axes2, axes3): 

149 axes.set_ylim(-0.5, len(targets_page)-0.5) 

150 axes.invert_yaxis() 

151 labelpos(axes, 2.5, 2.0) 

152 

153 axes1.axvline(0.0, color='black') 

154 

155 labels = [] 

156 lastpath = None 

157 for ipos, t in enumerate(targets_page): 

158 scodes = '.'.join(t.codes) 

159 if lastpath is None or lastpath != t.path: 

160 labels.append(t.path + '.' + scodes) 

161 if lastpath is not None: 

162 for axes in (axes1, axes2, axes3): 

163 axes.axhline(ipos-0.5, color='black') 

164 else: 

165 labels.append(scodes) 

166 

167 lastpath = t.path 

168 

169 for ipos, itarget in enumerate(range(ilo, ihi)): 

170 axes1.plot( 

171 residuals[-1, itarget], 

172 ipos, 

173 '|', 

174 ms=self.font_size, 

175 color='black') 

176 

177 if show_mean_residuals: 

178 axes1.plot( 

179 mean_residuals[itarget], 

180 ipos, 

181 's', 

182 ms=self.font_size, 

183 color='none', 

184 mec='black') 

185 

186 axes1.scatter( 

187 residuals[::10, itarget], 

188 util.num_full(icolor[::10].size, ipos, dtype=float), 

189 c=icolor[::10], 

190 cmap=cmap, 

191 norm=norm, 

192 alpha=0.5) 

193 

194 axes2.scatter( 

195 distances[::10, itarget]/km, 

196 util.num_full(icolor[::10].size, ipos, dtype=float), 

197 c=icolor[::10], 

198 cmap=cmap, 

199 norm=norm, 

200 alpha=0.5) 

201 

202 axes3.scatter( 

203 azimuths[::10, itarget], 

204 util.num_full(icolor[::10].size, ipos, dtype=float), 

205 c=icolor[::10], 

206 cmap=cmap, 

207 norm=norm, 

208 alpha=0.5) 

209 

210 axes1.set_yticks(num.arange(len(labels))) 

211 axes1.set_yticklabels(labels) 

212 axes1.set_xlabel('$T_{obs} - T_{syn}$ [s]') 

213 axes1.set_xlim(-residual_amax, residual_amax) 

214 

215 axes2.set_xlabel('Distance [km]') 

216 axes2.set_xlim(distance_min/km, distance_max/km) 

217 

218 axes3.set_xlabel('Azimuth [deg]') 

219 axes3.set_xlim(-180., 180.) 

220 

221 axes2.get_yaxis().set_ticks([]) 

222 axes2.get_yaxis().set_ticks([]) 

223 

224 axes3.get_yaxis().set_ticks([]) 

225 axes3.get_yaxis().set_ticks([]) 

226 

227 yield item, fig 

228 

229 

230class PickResidualsStationDistribution(StationDistributionPlot): 

231 ''' 

232 Plot showing residuals of seismic phase arrivals by azimuth and distance. 

233 ''' 

234 

235 name = 'pick_residuals_stations' 

236 

237 def make(self, environ): 

238 environ.setup_modelling() 

239 

240 cm = environ.get_plot_collection_manager() 

241 mpl_init(fontsize=self.font_size) 

242 

243 problem = environ.get_problem() 

244 dataset = environ.get_dataset() 

245 history = environ.get_history(subset='harvest') 

246 

247 cm.create_group_mpl( 

248 self, 

249 self.draw_figures(problem, dataset, history), 

250 title=u'Pick Residuals by Location', 

251 section='fits', 

252 feather_icon='target', 

253 description=u''' 

254Plot showing residuals for seismic phase arrivals by azimuth and distance. 

255 

256Station locations in dependence of distance and azimuth are shown. The center 

257of the plot corresponds to the origin of the search space. Optimized source 

258locations are shown as small black dots (subset of the solution ensemble). 

259''') 

260 

261 def draw_figures(self, problem, dataset, history): 

262 

263 target_index = {} 

264 i = 0 

265 for target in problem.targets: 

266 target_index[target] = i, i+target.nmisfits 

267 i += target.nmisfits 

268 

269 misfits = history.misfits[history.get_sorted_misfits_idx(), ...] 

270 gcms = problem.combine_misfits( 

271 misfits[:1, :, :], get_contributions=True)[0, :] 

272 models = history.get_sorted_primary_models()[::-1] 

273 

274 origin = problem.base_source 

275 

276 targets = [ 

277 t for t in problem.targets if isinstance(t, PhasePickTarget)] 

278 

279 ntargets = len(targets) 

280 nmodels = history.nmodels 

281 tts = num.zeros((nmodels, ntargets, 2)) 

282 

283 sdata = [] 

284 for imodel in range(nmodels): 

285 model = models[imodel, :] 

286 source = problem.get_source(model) 

287 sdata.append( 

288 (origin.distance_to(source), origin.azibazi_to(source)[0])) 

289 

290 results = problem.evaluate(model, targets=targets) 

291 for itarget, result in enumerate(results): 

292 result = results[itarget] 

293 tts[imodel, itarget, :] = result.tobs, result.tsyn 

294 

295 ok = num.all(num.isfinite(tts), axis=2) 

296 ok = num.all(ok, axis=0) 

297 

298 targets_ok = [ 

299 target for (itarget, target) in enumerate(targets) if ok[itarget]] 

300 tts = tts[:, ok, :] 

301 residuals = tts[:, :, 0] - tts[:, :, 1] 

302 mean_residuals = num.mean(residuals, axis=0) 

303 

304 rlo, rhi = num.percentile(mean_residuals, [10., 90.]) 

305 residual_amax = max(abs(rlo), abs(rhi)) 

306 

307 target_to_residual = dict( 

308 (target, residual) 

309 for (target, residual) 

310 in zip(targets_ok, mean_residuals)) 

311 

312 cg_to_targets = meta.gather(targets, lambda t: (t.path, )) 

313 

314 cgs = sorted(cg_to_targets.keys()) 

315 

316 for cg in cgs: 

317 cg_str = '.'.join(cg) 

318 

319 targets = cg_to_targets[cg] 

320 if len(targets) == 0: 

321 continue 

322 

323 assert all(target_index[target][0] == target_index[target][1] - 1 

324 for target in targets) 

325 

326 itargets = num.array( 

327 [target_index[target][0] for target in targets]) 

328 

329 labels = ['.'.join(x for x in t.codes[:3] if x) for t in targets] 

330 

331 azimuths = num.array( 

332 [origin.azibazi_to(t)[0] for t in targets]) 

333 distances = num.array( 

334 [t.distance_to(origin) for t in targets]) 

335 residuals = num.array( 

336 [target_to_residual.get(t, num.nan) for t in targets]) 

337 

338 item = PlotItem( 

339 name='picks_contributions_%s' % cg_str, 

340 title=u'Pick residuals and contributions (%s)' % cg_str, 

341 description=u'\n\nMarkers are scaled according to their ' 

342 u'misfit contribution for the globally best ' 

343 u'source model.') 

344 

345 fig, axes, legend = self.plot_station_distribution( 

346 azimuths, distances, 

347 gcms[itargets], 

348 labels, 

349 colors=residuals, 

350 cnorm=(-residual_amax, residual_amax), 

351 cmap='RdYlBu', 

352 clabel='$T_{obs} - T_{syn}$ [s]', 

353 scatter_kwargs=dict(alpha=1.0), 

354 legend_title='Contribution') 

355 

356 sources = [ 

357 problem.get_source(x) for x in models[::10]] 

358 

359 azimuths = num.array( 

360 [origin.azibazi_to(s)[0] for s in sources]) 

361 distances = num.array( 

362 [s.distance_to(origin) for s in sources]) 

363 

364 axes.plot( 

365 azimuths*d2r, distances, 'o', ms=2.0, color='black', alpha=0.3) 

366 

367 yield (item, fig) 

368 

369 

370def get_plot_classes(): 

371 return [ 

372 PickResidualsPlot, 

373 PickResidualsStationDistribution]