1import math
2import logging
4import numpy as num
6from matplotlib import pyplot as plt
8from pyrocko.guts import Float, Bool, Tuple
9from pyrocko import gf
10from pyrocko.plot import automap, mpl_init, beachball, mpl_color
12from grond import stats
13from grond.plot.collection import PlotItem
14from grond.plot.config import PlotConfig
16logger = logging.getLogger('grond.problem.double_sf.plot')
18guts_prefix = 'grond'
20km = 1e3
23class SFForcePlot(PlotConfig):
24 ''' Maps showing horizontal and vertical force
25 of the best double single force model '''
27 name = 'forces_double_singleforce'
29 size_cm = Tuple.T(
30 2, Float.T(),
31 default=(15., 15.),
32 help='width and length of the figure in cm')
33 show_topo = Bool.T(
34 default=False,
35 help='show topography')
36 show_grid = Bool.T(
37 default=True,
38 help='show the lat/lon grid')
39 show_rivers = Bool.T(
40 default=True,
41 help='show rivers on the map')
42 radius = Float.T(
43 optional=True,
44 help='radius of the map around campaign center lat/lon')
46 def make(self, environ):
47 cm = environ.get_plot_collection_manager()
48 history = environ.get_history(subset='harvest')
49 optimiser = environ.get_optimiser()
50 ds = environ.get_dataset()
52 environ.setup_modelling()
54 cm.create_group_automap(
55 self,
56 self.draw_best_sf(ds, history, optimiser),
57 title=u'Single Force Source Forces',
58 section='solution',
59 feather_icon='map',
60 description=u'''
61Maps show located force vectors of the best double Single Force Source model.
63Arrows show the modelled forces (red arrows). The top plot shows the horizontal
64forces and the bottom plot the vertical force. The dot indicates the location
65of the best double single force source model.
66''')
68 def draw_best_sf(self, ds, history, optimiser, vertical=False):
69 from grond.core import make_stats
71 source = history.get_best_source()
73 problem = history.problem
74 models = history.models
76 stats = make_stats(
77 problem, models, history.get_primary_chain_misfits())
79 def plot_double_sf(source, stats, ifig, vertical=False):
80 orient = 'vertical' if vertical else 'horizontal'
82 item = PlotItem(
83 name='fig_%i' % ifig,
84 attributes={},
85 title=u'Best %s double single force model force vector' % (
86 orient),
87 description=u'''
88Double single force source %s force vector for the best model (red). The circle
89shows the 95%% confidence ellipse. Orange points indicate the location of each
90individual single force, while the square shows the combined centroid location.
91''' % orient)
93 event = source.pyrocko_event()
95 radius = self.radius
96 if radius is None or radius < 30.*km:
97 logger.warn(
98 'Radius too small, defaulting to 30 km')
99 radius = 30*km
101 m = automap.Map(
102 width=self.size_cm[0],
103 height=self.size_cm[1],
104 lat=event.effective_lat,
105 lon=event.effective_lon,
106 radius=radius,
107 show_topo=self.show_topo,
108 show_grid=self.show_grid,
109 show_rivers=self.show_rivers,
110 color_wet=(216, 242, 254),
111 color_dry=(238, 236, 230))
113 source.lat, source.lon = event.effective_lat, event.effective_lon
115 if isinstance(source, gf.DoubleSFSource):
116 sf1, sf2 = source.split()
117 offset_scale = source.force
118 f = 'r'
119 elif isinstance(source, gf.CombiSFSource):
120 sf1, sf2 = source.subsources
122 sf1.lat, sf1.lon = source.lat, source.lon
123 sf2.lat, sf2.lon = source.lat, source.lon
125 offset_scale = num.max((sf1.force, sf2.force))
126 f = ''
128 size = num.linalg.norm(self.size_cm)
130 scale = (size / 5.) / offset_scale
132 stats_dict = stats.get_values_dict()
134 if vertical:
135 rows = [[
136 sf1.effective_lon, sf1.effective_lat,
137 0., -sf1.fd * scale,
138 (stats_dict[f + 'fn1.std'] + stats_dict[f + 'fe1.std']),
139 stats_dict[f + 'fd1.std'],
140 0.]]
142 rows.append([
143 sf2.effective_lon, sf2.effective_lat,
144 0., -sf2.fd * scale,
145 (stats_dict[f + 'fn2.std'] + stats_dict[f + 'fe2.std']),
146 stats_dict[f + 'fd2.std'],
147 0.])
149 else:
150 rows = [[
151 sf1.effective_lon, sf1.effective_lat,
152 sf1.fe * scale, sf1.fn * scale,
153 stats_dict[f + 'fe1.std'],
154 stats_dict[f + 'fn1.std'],
155 0.]]
157 rows.append([
158 sf2.effective_lon, sf2.effective_lat,
159 sf2.fe * scale, sf2.fn * scale,
160 stats_dict[f + 'fe2.std'],
161 stats_dict[f + 'fn2.std'],
162 0.])
164 fontsize = 10.
166 default_psxy_style = {
167 'h': 0,
168 'W': '2.0p,red',
169 'A': '+p4p,black+e+a40',
170 'G': 'red',
171 't': 30,
172 'L': True,
173 'S': 'e1c/0.95/%d' % fontsize,
174 }
176 m.gmt.psvelo(
177 in_rows=rows,
178 *m.jxyr,
179 **default_psxy_style)
181 m.gmt.psxy(
182 S='c8p',
183 in_rows=[
184 [sf.effective_lon, sf.effective_lat] for sf in (sf1, sf2)],
185 W='1p,black',
186 G='orange3',
187 *m.jxyr)
189 m.gmt.psxy(
190 S='s8p',
191 in_rows=[[source.effective_lon, source.effective_lat]],
192 W='1p,black',
193 G='slategray3',
194 *m.jxyr)
196 return (item, m)
198 ifig = 0
199 for vertical in (False, True):
200 yield plot_double_sf(source, stats, ifig, vertical)
201 ifig += 1
204class DoubleSFDecompositionPlot(PlotConfig):
205 '''
206 Double Single Force decomposition plot.
207 '''
209 name = 'sf_decomposition'
210 size_cm = Tuple.T(2, Float.T(), default=(15., 5.))
212 def make(self, environ):
213 cm = environ.get_plot_collection_manager()
214 history = environ.get_history(subset='harvest')
215 mpl_init(fontsize=self.font_size)
216 cm.create_group_mpl(
217 self,
218 self.draw_figures(history),
219 title=u'Single Force Decomposition',
220 section='solution',
221 feather_icon='sun',
222 description=u'''
223Double Single Force decomposition of the best-fitting solution into its single
224force components.
226Shown are the ensemble best and the ensemble mean. The symbol size indicates
227the relative strength of the components. The inversion result is consistent
228and stable if ensemble mean and ensemble
229best have similar symbol size and patterns.
230''')
232 def draw_figures(self, history):
234 fontsize = self.font_size
236 fig = plt.figure(figsize=self.size_inch)
237 axes = fig.add_subplot(1, 1, 1, aspect=1.0)
238 fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
240 problem = history.problem
241 models = history.models
243 if models.size == 0:
244 logger.warn('Empty models vector.')
245 return []
247 # ref_source = problem.base_source
249 mean_source = stats.get_mean_source(
250 problem, history.models)
252 best_source = history.get_best_source()
254 nlines_max = int(round(self.size_cm[1] / 5. * 4. - 1.0))
256 def get_deco(source):
257 if isinstance(source, gf.DoubleSFSource):
258 return [source] + source.split()
259 elif isinstance(source, gf.CombiSFSource):
260 sf1, sf2 = source.subsources
261 f1 = num.array([sf1.fn, sf1.fe, sf1.fd])
262 f2 = num.array([sf2.fn, sf2.fe, sf2.fd])
264 force = num.linalg.norm([f1 + f2])
266 mix = sf2.force / force
268 source = gf.DoubleSFSource(
269 force=force,
270 rfn1=sf1.fn / force,
271 rfe1=sf1.fe / force,
272 rfd1=sf1.fd / force,
273 rfn2=sf2.fn / force,
274 rfe2=sf2.fe / force,
275 rfd2=sf2.fd / force,
276 mix=mix)
278 return [source, sf1, sf2]
280 lines = []
281 lines.append(
282 ('Ensemble best', get_deco(best_source), mpl_color('aluminium5')))
284 lines.append(
285 ('Ensemble mean', get_deco(mean_source), mpl_color('aluminium5')))
287 force_max = max(sf.force for (_, line, _) in lines for sf in line)
289 for xpos, label in [
290 (0., 'Double SF'),
291 (2., 'SF 1'),
292 (4., 'SF 2')]:
294 axes.annotate(
295 label,
296 xy=(1 + xpos, nlines_max),
297 xycoords='data',
298 xytext=(0., 0.),
299 textcoords='offset points',
300 ha='center',
301 va='center',
302 color='black',
303 fontsize=fontsize)
305 for i, (label, deco, color_t) in enumerate(lines):
306 ypos = nlines_max - i - 1.0
308 [dsf, sf1, sf2] = deco
310 size0 = dsf.force / force_max
312 axes.annotate(
313 label,
314 xy=(-2., ypos),
315 xycoords='data',
316 xytext=(0., 0.),
317 textcoords='offset points',
318 ha='left',
319 va='center',
320 color='black',
321 fontsize=fontsize)
323 for xpos, sf_part, ratio, ops in [
324 (0., dsf, 1., '='),
325 (2., sf1, sf1.force / force_max, '+'),
326 (4., sf2, sf2.force / force_max, None)]:
328 if ratio > 1e-4:
329 try:
330 beachball.plot_singleforce_beachball_mpl(
331 sf_part.fn, sf_part.fe, sf_part.fd, axes,
332 position=(1. + xpos, ypos),
333 size=0.9 * size0 * math.sqrt(ratio),
334 size_units='data',
335 color_t=color_t,
336 linewidth=1.0)
338 except beachball.BeachballError as e:
339 logger.warn(str(e))
341 axes.annotate(
342 'ERROR',
343 xy=(1. + xpos, ypos),
344 ha='center',
345 va='center',
346 color='red',
347 fontsize=fontsize)
349 else:
350 axes.annotate(
351 'N/A',
352 xy=(1. + xpos, ypos),
353 ha='center',
354 va='center',
355 color='black',
356 fontsize=fontsize)
358 if ops is not None:
359 axes.annotate(
360 ops,
361 xy=(2. + xpos, ypos),
362 ha='center',
363 va='center',
364 color='black',
365 fontsize=fontsize)
367 axes.axison = False
368 axes.set_xlim(-2.25, 9.75)
369 axes.set_ylim(-0.5, nlines_max+0.5)
371 item = PlotItem(name='main')
372 return [[item, fig]]