Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/dynamic_rupture/plot.py: 23%

1import numpy as num 

2import matplotlib.pyplot as plt 

3from matplotlib.ticker import FuncFormatter 

4from matplotlib import colors, patheffects 

5 

6from pyrocko.guts import Tuple, Float, Bool 

7from pyrocko.gf.seismosizer import map_anchor 

8from pyrocko.plot import mpl_init, mpl_graph_color, mpl_color 

9from pyrocko.plot.dynamic_rupture import RuptureMap 

10 

11from grond.plot.config import PlotConfig 

12from grond.plot.collection import PlotItem 

13 

14km = 1e3 

15 

16 

17def km_fmt(x, p): 

18 return '%.1f' % (x / km) 

19 

20 

21def get_source(history, source_type='mean'): 

22 if source_type == 'mean': 

23 return history.get_mean_source() 

24 elif source_type == 'best': 

25 return history.get_best_source() 

26 else: 

27 raise ValueError('current source_type is not defined.') 

28 

29 

30class DynamicRuptureSlipMap(PlotConfig): 

31 ''' 

32 Slip map of best solution. 

33 ''' 

34 name = 'rupture_slip_map' 

35 dt_contour = Float.T( 

36 default=0.5, 

37 help='Rupture propagation contourline interval in seconds.') 

38 size_cm = Tuple.T(2, Float.T(), default=(20., 20.)) 

39 

40 def make(self, environ): 

41 environ.setup_modelling() 

42 cm = environ.get_plot_collection_manager() 

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

44 problem = environ.get_problem() 

45 

46 mpl_init(fontsize=self.font_size) 

47 cm.create_group_mpl( 

48 self, 

49 self.draw_figures(history, problem), 

50 title=u'Slip Map', 

51 section='solution', 

52 feather_icon='grid', # alternatively: wind 

53 description=u''' 

54Slip distribution and rake of the pseudo dynamic rupture model. The absolute 

55slip is color coded. The arrows are scaled based on the in plane slip. The 

56black star marks the nucleation point. Contour lines indicate the rupture 

57evolution in %.1f s intervals. 

58''' % self.dt_contour) 

59 

60 def draw_figures(self, history, problem): 

61 store_ids = problem.get_gf_store_ids() 

62 store = problem.get_gf_store(store_ids[0]) 

63 

64 interpolation = 'nearest_neighbor' 

65 

66 for i, plabel in enumerate(('best', 'mean')): 

67 source = get_source(history, source_type=plabel) 

68 anchor_x, anchor_y = map_anchor[source.anchor] 

69 

70 fig, ax = plt.subplots(1, 1) 

71 

72 # ToDo in function with "mean", "best" as arg 

73 source.discretize_patches(store, interpolation) 

74 patches = source.patches 

75 dislocations = source.get_slip(scale_slip=True) 

76 

77 patches_x = num.array([p.ix for p in patches])\ 

78 .reshape(source.nx, source.ny) 

79 patches_y = num.array([p.iy for p in patches])\ 

80 .reshape(source.nx, source.ny) 

81 patches_t = num.array([p.time for p in patches])\ 

82 .reshape(source.nx, source.ny) 

83 

84 patches_x += (anchor_x + 1.) / 2. * source.length 

85 patches_y += (anchor_y + 1.) / 2. * source.width 

86 

87 abs_disloc = num.linalg.norm(dislocations, axis=1) 

88 abs_disloc = abs_disloc.reshape(source.nx, source.ny) 

89 

90 im = ax.imshow( 

91 abs_disloc.T, 

92 cmap='YlOrRd', 

93 origin='upper', 

94 aspect='equal', 

95 extent=(0., source.length, 0., source.width)) 

96 

97 patches_t -= patches_t.min() 

98 

99 nlevels = patches_t.max() // self.dt_contour 

100 contours = num.arange(nlevels) * self.dt_contour 

101 contours += patches_t.min() 

102 

103 def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100): 

104 return colors.LinearSegmentedColormap.from_list( 

105 'trunc({n},{a:.2f},{b:.2f})'.format( 

106 n=cmap.name, a=minval, b=maxval), 

107 cmap(num.linspace(minval, maxval, n))) 

108 

109 cmap = truncate_colormap(plt.get_cmap('winter'), 0., 0.8) 

110 

111 contour = ax.contour( 

112 patches_x, patches_y, patches_t, 

113 levels=contours, alpha=.8, colors='k') 

114 

115 labels = ax.clabel( 

116 contour, contour.levels[::2], 

117 inline=True, fmt='%.1f s') 

118 

119 for label in labels: 

120 label.set_rotation(0.) 

121 label.set_fontweight('semibold') 

122 label.set_fontsize('small') 

123 label.set_path_effects([ 

124 patheffects.Stroke(linewidth=1.25, foreground='beige'), 

125 patheffects.Normal()]) 

126 

127 slip_strike = dislocations[:, 0] 

128 slip_dip = dislocations[:, 1] 

129 

130 patch_length = source.length / source.nx 

131 patch_width = source.width / source.ny 

132 

133 x_quiver = num.repeat(num.arange(source.nx), source.ny)\ 

134 * patch_length + patch_length/2. 

135 y_quiver = num.tile(num.arange(source.ny), source.nx)\ 

136 * patch_width + patch_width/2. 

137 

138 ax.quiver( 

139 x_quiver, y_quiver, slip_strike, slip_dip, 

140 facecolor='none', edgecolor='k', linewidth=.7, 

141 scale=75 / abs_disloc.max(), headwidth=3, 

142 cmap='YlOrRd', alpha=.6) 

143 

144 ax.invert_yaxis() 

145 

146 nucleation_x = ((source.nucleation_x + 1.) / 2.) * source.length 

147 nucleation_y = ((source.nucleation_y + 1.) / 2.) * source.width 

148 

149 ax.scatter( 

150 nucleation_x, nucleation_y, 

151 s=60, color='w', edgecolor='k', marker='*', 

152 linewidths=.5, alpha=.7) 

153 

154 ax.xaxis.set_major_formatter(FuncFormatter(km_fmt)) 

155 ax.yaxis.set_major_formatter(FuncFormatter(km_fmt)) 

156 

157 ax.set_xlabel('Length [km]') 

158 ax.set_ylabel('Width [km]') 

159 

160 cmap = fig.colorbar( 

161 im, orientation='horizontal', pad=0.2, shrink=.8, 

162 ax=ax, format='%.2f') 

163 cmap.set_label('Slip [m]') 

164 

165 item = PlotItem( 

166 name='_'.join(['ensemble', plabel.lower()]), 

167 title='ensemble %s slip distribution' % (plabel.lower()), 

168 description=u'') 

169 

170 yield item, fig 

171 

172 

173class DynamicRuptureSTF(PlotConfig): 

174 ''' 

175 Slip map of best solution. 

176 ''' 

177 name = 'rupture_source_time_function' 

178 dt_sampling = Float.T( 

179 optional=True, 

180 help='Moment rate sampling interval in seconds') 

181 size_cm = Tuple.T(2, Float.T(), default=(15., 15.)) 

182 

183 def _get_store(self, problem): 

184 store_ids = problem.get_gf_store_ids() 

185 

186 stores = [problem.get_gf_store(store_id) for store_id in store_ids] 

187 dt = num.array([store.config.deltat for store in stores]) 

188 min_index = num.where(dt == min(dt))[0] 

189 

190 return stores[min_index[0]] 

191 

192 def make(self, environ): 

193 environ.setup_modelling() 

194 cm = environ.get_plot_collection_manager() 

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

196 problem = environ.get_problem() 

197 

198 store = self._get_store(problem=problem) 

199 

200 dt_sampling = store.config.deltat if self.dt_sampling is None \ 

201 else self.dt_sampling 

202 

203 mpl_init(fontsize=self.font_size) 

204 cm.create_group_mpl( 

205 self, 

206 self.draw_figures(history, problem), 

207 title=u'Moment Rate Function (STF)', 

208 section='solution', 

209 feather_icon='zap', # alternatively: activity 

210 description=u''' 

211Source time function (moment release) of the pseudo dynamic rupture model. 

212The moment rate function is sampled in %.1f s intervals. 

213''' % dt_sampling) 

214 

215 def draw_figures(self, history, problem): 

216 store = self._get_store(problem=problem) 

217 

218 for i, plabel in enumerate(('best', 'mean')): 

219 source = get_source(history, source_type=plabel) 

220 

221 fig, ax = plt.subplots(1, 1) 

222 

223 dt_sampling = store.config.deltat if self.dt_sampling is None \ 

224 else self.dt_sampling 

225 

226 mrate, times = source.get_moment_rate( 

227 store=store, deltat=dt_sampling) 

228 

229 try: 

230 mrate_max = mrate.max() 

231 

232 ax.fill_between( 

233 times, 

234 mrate / mrate_max, 

235 color=mpl_color(x='aluminium2'), 

236 alpha=0.7) 

237 

238 ax.scatter( 

239 times, 

240 mrate / mrate_max, 

241 c=mpl_graph_color(0)) 

242 

243 ax.set_xlabel('Time [s]') 

244 ax.set_ylabel(r'$\dot{M}$ / %.2e Nm/s' % mrate_max) 

245 

246 ax.set_xlim([0, times.max() + dt_sampling]) 

247 ax.grid(True) 

248 

249 except ValueError: 

250 pass 

251 

252 item = PlotItem( 

253 name='_'.join(['ensemble', plabel.lower()]), 

254 title='ensemble %s source time function' % (plabel.lower()), 

255 description=u'') 

256 

257 yield item, fig 

258 

259 

260class DynamicRuptureMap(PlotConfig): 

261 ''' 

262 Overview map of best solution. 

263 ''' 

264 name = 'rupture_overview_map' 

265 size_cm = Tuple.T( 

266 2, Float.T(), 

267 default=(20., 20.), 

268 help='width and length of the figure in cm') 

269 show_topo = Bool.T( 

270 default=False, 

271 help='show topography') 

272 show_grid = Bool.T( 

273 default=True, 

274 help='show the lat/lon grid') 

275 show_rivers = Bool.T( 

276 default=True, 

277 help='show rivers on the map') 

278 radius = Float.T( 

279 optional=True, 

280 help='radius of the map around campaign center lat/lon') 

281 

282 def make(self, environ): 

283 environ.setup_modelling() 

284 cm = environ.get_plot_collection_manager() 

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

286 problem = environ.get_problem() 

287 

288 cm.create_group_automap( 

289 self, 

290 self.draw_rupture_map(history, problem), 

291 title=u'Rupture Dislocations', 

292 section='solution', 

293 feather_icon='map', 

294 description=u''' 

295Maps showing orientation and dislocation of the PseudoDynamicRupture. 

296''') 

297 

298 def draw_rupture_map(self, history, problem): 

299 from pyrocko import orthodrome as pod 

300 

301 store_ids = problem.get_gf_store_ids() 

302 store = problem.get_gf_store(store_ids[0]) 

303 

304 def plot_rupture(source, model_name): 

305 item = PlotItem( 

306 name='ensemble_%s' % model_name, 

307 attributes={}, 

308 title=u'Final static rupture dislocation - %s model' 

309 % model_name, 

310 description=u''' 

311Static rupture dislocation from %s model.''' % model_name) 

312 

313 lat, lon = pod.ne_to_latlon( 

314 source.lat, 

315 source.lon, 

316 source.north_shift, 

317 source.east_shift) 

318 

319 map_kwargs = dict( 

320 lat=lat, 

321 lon=lon, 

322 radius=self.radius or source.length, 

323 width=self.size_cm[0], 

324 height=self.size_cm[1], 

325 source=source, 

326 show_topo=self.show_topo, 

327 show_grid=self.show_grid, 

328 show_rivers=self.show_rivers, 

329 color_wet=(216, 242, 254), 

330 color_dry=(238, 236, 230)) 

331 

332 source.discretize_patches(store, interpolation='nearest_neighbor') 

333 

334 cmap = 'summer' 

335 

336 m = RuptureMap(**map_kwargs) 

337 m.draw_dislocation(cmap=cmap) 

338 m.draw_time_contour(store) 

339 m.draw_nucleation_point() 

340 m.draw_top_edge() 

341 

342 return (item, m) 

343 

344 for i, plabel in enumerate(('best', 'mean')): 

345 source = get_source(history, source_type=plabel) 

346 

347 yield plot_rupture(source, plabel)