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 19:53 +0000
1import math
3from matplotlib import pyplot as plt, colors as mcolors
4import numpy as num
6from pyrocko import util
7from pyrocko.guts import Tuple, Float
8from pyrocko.plot import mpl_init, mpl_margins
10from grond import meta
11from grond.plot.config import PlotConfig
12from grond.plot.collection import PlotItem
14from .target import PhasePickTarget
15from ..plot import StationDistributionPlot
17km = 1000.
18d2r = math.pi / 180.
21class PickResidualsPlot(PlotConfig):
22 '''Travel time residuals plot.'''
24 name = 'pick_residuals'
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]')
31 def make(self, environ):
32 environ.setup_modelling()
34 cm = environ.get_plot_collection_manager()
35 mpl_init(fontsize=self.font_size)
37 ds = environ.get_dataset()
38 history = environ.get_history(subset='harvest')
39 optimiser = environ.get_optimiser()
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.
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''')
56 def draw_figures(self, ds, history, optimiser):
57 show_mean_residuals = False
59 # gms = history.get_sorted_primary_misfits()[::-1]
60 models = history.get_sorted_primary_models()[::-1]
61 problem = history.problem
63 targets = [
64 t for t in problem.targets if
65 isinstance(t, PhasePickTarget)]
67 # targets.sort(key=lambda t: t.distance_to(problem.base_source))
68 # tpaths = sorted(set(t.path for t in targets))
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))
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])
86 azimuths[imodel, itarget] = source.azibazi_to(
87 targets[itarget])[0]
89 ok = num.all(num.isfinite(tts), axis=2)
90 ok = num.all(ok, axis=0)
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)
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)
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]
116 icolor = num.arange(nmodels)
118 norm = mcolors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor))
119 cmap = plt.get_cmap('coolwarm')
121 napprox = max(
122 1, int(round((self.size_points[1] / self.font_size - 7))))
124 npages = (ntargets-1) // napprox + 1
125 ntargets_page = (ntargets-1) // npages + 1
127 for ipage in range(npages):
129 ilo, ihi = ipage*ntargets_page, (ipage+1)*ntargets_page
130 ihi = min(len(targets), ihi)
132 targets_page = targets[ilo:ihi]
133 item = PlotItem(
134 name='fig_%02i' % ipage)
136 item.attributes['targets'] = [t.string_id() for t in targets_page]
138 fig = plt.figure(figsize=self.size_inch)
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)
144 axes1 = fig.add_subplot(1, 3, 1)
145 axes2 = fig.add_subplot(1, 3, 2)
146 axes3 = fig.add_subplot(1, 3, 3)
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)
153 axes1.axvline(0.0, color='black')
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)
167 lastpath = t.path
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')
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')
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)
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)
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)
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)
215 axes2.set_xlabel('Distance [km]')
216 axes2.set_xlim(distance_min/km, distance_max/km)
218 axes3.set_xlabel('Azimuth [deg]')
219 axes3.set_xlim(-180., 180.)
221 axes2.get_yaxis().set_ticks([])
222 axes2.get_yaxis().set_ticks([])
224 axes3.get_yaxis().set_ticks([])
225 axes3.get_yaxis().set_ticks([])
227 yield item, fig
230class PickResidualsStationDistribution(StationDistributionPlot):
231 '''
232 Plot showing residuals of seismic phase arrivals by azimuth and distance.
233 '''
235 name = 'pick_residuals_stations'
237 def make(self, environ):
238 environ.setup_modelling()
240 cm = environ.get_plot_collection_manager()
241 mpl_init(fontsize=self.font_size)
243 problem = environ.get_problem()
244 dataset = environ.get_dataset()
245 history = environ.get_history(subset='harvest')
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.
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''')
261 def draw_figures(self, problem, dataset, history):
263 target_index = {}
264 i = 0
265 for target in problem.targets:
266 target_index[target] = i, i+target.nmisfits
267 i += target.nmisfits
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]
274 origin = problem.base_source
276 targets = [
277 t for t in problem.targets if isinstance(t, PhasePickTarget)]
279 ntargets = len(targets)
280 nmodels = history.nmodels
281 tts = num.zeros((nmodels, ntargets, 2))
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]))
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
295 ok = num.all(num.isfinite(tts), axis=2)
296 ok = num.all(ok, axis=0)
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)
304 rlo, rhi = num.percentile(mean_residuals, [10., 90.])
305 residual_amax = max(abs(rlo), abs(rhi))
307 target_to_residual = dict(
308 (target, residual)
309 for (target, residual)
310 in zip(targets_ok, mean_residuals))
312 cg_to_targets = meta.gather(targets, lambda t: (t.path, ))
314 cgs = sorted(cg_to_targets.keys())
316 for cg in cgs:
317 cg_str = '.'.join(cg)
319 targets = cg_to_targets[cg]
320 if len(targets) == 0:
321 continue
323 assert all(target_index[target][0] == target_index[target][1] - 1
324 for target in targets)
326 itargets = num.array(
327 [target_index[target][0] for target in targets])
329 labels = ['.'.join(x for x in t.codes[:3] if x) for t in targets]
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])
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.')
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')
356 sources = [
357 problem.get_source(x) for x in models[::10]]
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])
364 axes.plot(
365 azimuths*d2r, distances, 'o', ms=2.0, color='black', alpha=0.3)
367 yield (item, fig)
370def get_plot_classes():
371 return [
372 PickResidualsPlot,
373 PickResidualsStationDistribution]