Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/singleforce/plot.py: 23%
125 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-06-12 14:01 +0000
1import logging
3import numpy as num
4from matplotlib import cm, patches
6from pyrocko import orthodrome as pod
7from pyrocko import gf
8from pyrocko.guts import Float, Bool, Tuple
10from pyrocko.plot import mpl_color, mpl_init, automap
12from grond.plot.section import SectionPlotConfig, SectionPlot
13from grond.plot.collection import PlotItem
14from grond.plot.config import PlotConfig
15from grond.problems.plot import fixlim
16from matplotlib import pyplot as plt
18logger = logging.getLogger('grond.problem.singleforce.plot')
20guts_prefix = 'grond'
22km = 1e3
25class SFLocationPlot(SectionPlotConfig):
26 ''' MT location plot of the best solutions in three cross-sections. '''
27 name = 'location_sf'
28 normalisation_gamma = Float.T(
29 default=3.,
30 help='Normalisation of colors and alpha as :math:`x^\\gamma`.'
31 'A linear colormap/alpha with :math:`\\gamma=1`.')
33 def make(self, environ):
34 environ.setup_modelling()
35 cm = environ.get_plot_collection_manager()
36 history = environ.get_history(subset='harvest')
37 mpl_init(fontsize=self.font_size)
38 self._to_be_closed = []
39 cm.create_group_mpl(
40 self,
41 self.draw_figures(history),
42 title=u'SF Location',
43 section='solution',
44 feather_icon='target',
45 description=u'''
46Location plot of the ensemble of best solutions in three cross-sections.
48The coordinate range is defined by the search space given in the config file.
49Dots show best single force mechanisms, and colors indicate low (red) and
50high (blue) misfit.
51''')
52 for obj in self._to_be_closed:
53 obj.close()
55 def draw_figures(self, history, color_p_axis=False):
56 from matplotlib import colors
58 color = 'black'
59 fontsize = self.font_size
60 markersize = fontsize * 1.5
62 problem = history.problem
63 sp = SectionPlot(config=self)
64 self._to_be_closed.append(sp)
66 fig = sp.fig
67 axes_en = sp.axes_xy
68 axes_dn = sp.axes_zy
69 axes_ed = sp.axes_xz
71 bounds = problem.get_combined_bounds()
73 models = history.get_sorted_primary_models()[::-1]
75 iorder = num.arange(history.nmodels)
77 for parname, set_label, set_lim in [
78 ['east_shift', sp.set_xlabel, sp.set_xlim],
79 ['north_shift', sp.set_ylabel, sp.set_ylim],
80 ['depth', sp.set_zlabel, sp.set_zlim]]:
82 ipar = problem.name_to_index(parname)
83 par = problem.combined[ipar]
84 set_label(par.get_label())
85 xmin, xmax = fixlim(*par.scaled(bounds[ipar]))
86 set_lim(xmin, xmax)
88 def scale_size(source):
89 return markersize * 1.5
91 for axes, xparname, yparname in [
92 (axes_en, 'east_shift', 'north_shift'),
93 (axes_dn, 'depth', 'north_shift'),
94 (axes_ed, 'east_shift', 'depth')]:
96 ixpar = problem.name_to_index(xparname)
97 iypar = problem.name_to_index(yparname)
99 xpar = problem.combined[ixpar]
100 ypar = problem.combined[iypar]
102 xmin, xmax = fixlim(*xpar.scaled(bounds[ixpar]))
103 ymin, ymax = fixlim(*ypar.scaled(bounds[iypar]))
105 try:
106 axes.set_facecolor(mpl_color('aluminium1'))
107 except AttributeError:
108 axes.patch.set_facecolor(mpl_color('aluminium1'))
110 rect = patches.Rectangle(
111 (xmin, ymin), xmax-xmin, ymax-ymin,
112 facecolor=mpl_color('white'),
113 edgecolor=mpl_color('aluminium2'),
114 zorder=1)
116 axes.add_patch(rect)
118 cmap = cm.ScalarMappable(
119 norm=colors.PowerNorm(
120 gamma=self.normalisation_gamma,
121 vmin=iorder.min(),
122 vmax=iorder.max()),
124 cmap=plt.get_cmap('coolwarm'))
126 for ix, x in enumerate(models):
128 source = problem.get_source(x)
129 fx = problem.extract(x, ixpar)
130 fy = problem.extract(x, iypar)
131 sx, sy = xpar.scaled(fx), ypar.scaled(fy)
133 # TODO: Add rotation in cross-sections
134 color = cmap.to_rgba(iorder[ix])
136 alpha = (iorder[ix] - iorder.min()) / \
137 float(iorder.max() - iorder.min())
138 alpha = alpha**self.normalisation_gamma
140 axes.scatter(
141 [sx], [sy],
142 c=[color],
143 s=[scale_size(source)],
144 alpha=alpha,
145 zorder=2)
147 item = PlotItem(name='main')
148 return [[item, fig]]
151class SFForcePlot(PlotConfig):
152 ''' Maps showing horizontal and vertical force
153 of the best single force model '''
155 name = 'forces_singleforce'
157 size_cm = Tuple.T(
158 2, Float.T(),
159 default=(15., 15.),
160 help='width and length of the figure in cm')
161 show_topo = Bool.T(
162 default=False,
163 help='show topography')
164 show_grid = Bool.T(
165 default=True,
166 help='show the lat/lon grid')
167 show_rivers = Bool.T(
168 default=True,
169 help='show rivers on the map')
170 radius = Float.T(
171 optional=True,
172 help='radius of the map around campaign center lat/lon')
174 def make(self, environ):
175 cm = environ.get_plot_collection_manager()
176 history = environ.get_history(subset='harvest')
177 optimiser = environ.get_optimiser()
178 ds = environ.get_dataset()
180 environ.setup_modelling()
182 cm.create_group_automap(
183 self,
184 self.draw_best_sf(ds, history, optimiser),
185 title=u'Single Force Source Forces',
186 section='solution',
187 feather_icon='map',
188 description=u'''
189Maps showing location and force vectors of the best Single Force Source model.
191Arrows show the modelled forces (red arrows). The top plot shows the horizontal
192forces and the bottom plot the vertical force. The dot indicates the location
193of the best single force source model.
194''')
196 def draw_best_sf(self, ds, history, optimiser, vertical=False):
197 from grond.core import make_stats
199 source = history.get_best_source()
201 problem = history.problem
202 models = history.models
204 stats = make_stats(
205 problem, models, history.get_primary_chain_misfits())
207 def plot_sf(source, stats, ifig, vertical=False):
208 orient = 'vertical' if vertical else 'horizontal'
210 item = PlotItem(
211 name='fig_%i' % ifig,
212 attributes={},
213 title=u'Best %s single force model force vector' % orient,
214 description=u'''
215Single force source %s force vector for the best model (red). The circle shows
216the 95%% confidence ellipse.
217''' % orient)
219 event = source.pyrocko_event()
221 radius = self.radius
222 if radius is None or radius < 100.*km:
223 logger.warn(
224 'Radius too small, defaulting to 100 km')
225 radius = 100*km
227 m = automap.Map(
228 width=self.size_cm[0],
229 height=self.size_cm[1],
230 lat=event.lat,
231 lon=event.lon,
232 radius=radius,
233 show_topo=self.show_topo,
234 show_grid=self.show_grid,
235 show_rivers=self.show_rivers,
236 color_wet=(216, 242, 254),
237 color_dry=(238, 236, 230))
239 if isinstance(source, gf.SFSource):
240 components = ['fn', 'fe', 'fd']
241 elif isinstance(source, gf.SimpleLandslideSource):
242 components = ['impulse_n', 'impulse_e', 'impulse_d']
244 offset_scale = num.abs(
245 [getattr(source, comp) for comp in components]).sum()
246 size = num.linalg.norm(self.size_cm)
248 scale = (size / 5.) / offset_scale
250 lat, lon = pod.ne_to_latlon(
251 event.lat,
252 event.lon,
253 source.north_shift,
254 source.east_shift)
256 stats_dict = stats.get_values_dict()
258 if vertical:
259 rows = [
260 lon, lat,
261 0., -getattr(source, components[2]) * scale,
262 num.sqrt(
263 stats_dict[components[0]+'.std']**2
264 + stats_dict[components[1]+'.std']**2) * scale,
265 stats_dict[components[2]+'.std'] * scale,
266 0.]
268 else:
269 rows = [
270 lon, lat,
271 getattr(source, components[1]) * scale,
272 getattr(source, components[0]) * scale,
273 stats_dict[components[1]+'.std'] * scale,
274 stats_dict[components[0]+'.std'] * scale,
275 0.]
277 fontsize = 10.
279 default_psxy_style = {
280 'h': 0,
281 'W': '2.0p,red',
282 'A': '+p4p,black+e+a40',
283 'G': 'red',
284 't': 30,
285 'L': True,
286 'S': 'e1c/0.95/%d' % fontsize,
287 }
289 m.gmt.psvelo(
290 in_rows=[rows],
291 *m.jxyr,
292 **default_psxy_style)
294 m.gmt.psxy(
295 S='c10p',
296 in_rows=[[lon, lat]],
297 W='1p,black',
298 G='orange3',
299 *m.jxyr)
301 return (item, m)
303 ifig = 0
304 for vertical in (False, True):
305 yield plot_sf(source, stats, ifig, vertical)
306 ifig += 1
308# TODO Fuzzy Single Force Plot