Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/directivity.py: 87%

208 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import logging 

7import numpy as num 

8import matplotlib.pyplot as plt 

9 

10from matplotlib.cm import ScalarMappable 

11from matplotlib.ticker import FuncFormatter 

12 

13from pyrocko.plot import beachball 

14from pyrocko.gf.meta import Timing 

15from pyrocko.gf import LocalEngine, Target, RectangularSource, map_anchor 

16 

17 

18km = 1e3 

19r2d = 180. / num.pi 

20d2r = num.pi / 180. 

21 

22logger = logging.getLogger(__name__) 

23 

24 

25QUANTITY_LABEL = { 

26 'displacement': 'Displacement [m]', 

27 'velocity': 'Velocity [m/s]', 

28 'acceleration': 'Acceleration [m/s²]' 

29} 

30 

31 

32def get_azimuthal_targets( 

33 store_id, source, radius, 

34 azi_begin=0., azi_end=360., dazi=1., 

35 depth=0.0, 

36 interpolation='multilinear', 

37 components='RTZ', quantity='displacement'): 

38 

39 assert dazi > 0. 

40 assert azi_begin < azi_end 

41 

42 nstations = int((azi_end - azi_begin) // dazi) 

43 assert nstations > 0 

44 

45 azimuths = num.linspace(azi_begin, azi_end, nstations) 

46 

47 coords = num.zeros((2, nstations)) 

48 coords[0, :] = num.cos(azimuths*d2r) 

49 coords[1, :] = num.sin(azimuths*d2r) 

50 coords *= radius 

51 

52 dips = {'R': 0., 'T': 0., 'Z': -90.} 

53 for comp in components: 

54 assert comp in dips.keys() 

55 

56 target_kwargs = dict( 

57 quantity=quantity, 

58 interpolation=interpolation, 

59 store_id=store_id) 

60 

61 targets = [ 

62 Target( 

63 lat=source.lat, 

64 lon=source.lon, 

65 north_shift=coords[0, iazi] + source.north_shift, 

66 east_shift=coords[1, iazi] + source.east_shift, 

67 depth=depth, 

68 azimuth={ 

69 'R': azi, 

70 'T': azi+90., 

71 'Z': 0. 

72 }[channel], 

73 dip=dips[channel], 

74 codes=('', 'S%01d' % iazi, '', channel), 

75 **target_kwargs) 

76 for iazi, azi in enumerate(azimuths) 

77 for channel in components] 

78 

79 for target, azi in zip(targets, azimuths): 

80 target.azimuth = azi 

81 target.dazi = dazi 

82 

83 return targets, azimuths 

84 

85 

86def get_seismogram_array( 

87 response, fmin=None, fmax=None, 

88 component='R', envelope=False): 

89 resp = response 

90 assert len(resp.request.sources) == 1, 'more than one source in response' 

91 

92 tmin = None 

93 tmax = None 

94 traces = [] 

95 

96 for _, target, tr in response.iter_results(): 

97 if target.codes[-1] != component: 

98 continue 

99 assert hasattr(target, 'azimuth') 

100 assert target.dazi 

101 

102 if fmin is not None: 

103 tr.highpass(4, fmin, demean=False) 

104 

105 if fmax is not None: 

106 tr.lowpass(4, fmax, demean=False) 

107 

108 tmin = min(tmin, tr.tmin) if tmin else tr.tmin 

109 tmax = max(tmax, tr.tmax) if tmax else tr.tmax 

110 traces.append(tr) 

111 

112 for tr in traces: 

113 tr.extend(tmin, tmax, fillmethod='repeat') 

114 if envelope: 

115 tr.envelope() 

116 

117 data = num.array([tr.get_ydata() for tr in traces]) 

118 nsamples = data.shape[1] 

119 return data, num.linspace(tmin, tmax, nsamples) 

120 

121 

122def hillshade(array, azimuth, angle_altitude): 

123 azimuth = 360.0 - azimuth 

124 azi = azimuth * r2d 

125 alt = angle_altitude * r2d 

126 

127 x, y = num.gradient(array) 

128 slope = num.pi/2. - num.arctan(num.sqrt(x*x + y*y)) 

129 aspect = num.arctan2(-x, y) 

130 

131 shaded = num.sin(alt)*num.sin(slope) \ 

132 + num.cos(alt)*num.cos(slope)*num.cos((azi - num.pi/2.) - aspect) 

133 

134 return (shaded + 1.)/2. 

135 

136 

137def hillshade_seismogram_array( 

138 seismogram_array, rgba_map, 

139 shad_lim=(.4, .98), contrast=1., blend_mode='multiply'): 

140 assert blend_mode in ('multiply', 'screen'), 'unknown blend mode' 

141 assert shad_lim[0] < shad_lim[1], 'bad shading limits' 

142 from scipy.ndimage import convolve as im_conv 

143 # Light source from somewhere above - psychologically the best choice 

144 # from upper left 

145 ramp = num.array([[1., 0.], [0., -1.]]) * contrast 

146 

147 # convolution of two 2D arrays 

148 shad = im_conv(seismogram_array, ramp.T).ravel() 

149 shad *= -1. 

150 

151 # if there are strong artifical edges in the data, shades get 

152 # dominated by them. Cutting off the largest and smallest 2% of 

153 # # shades helps 

154 percentile2 = num.percentile(shad, 2.0) 

155 percentile98 = num.percentile(shad, 98.0) 

156 

157 shad[shad > percentile98] = percentile98 

158 shad[shad < percentile2] = percentile2 

159 

160 # # normalize shading 

161 shad -= num.nanmin(shad) 

162 shad /= num.nanmax(shad) 

163 

164 # # reduce range to balance gray color 

165 shad *= shad_lim[1] - shad_lim[0] 

166 shad += shad_lim[0] 

167 

168 if blend_mode == 'screen': 

169 rgba_map[:, :3] = 1. - ((1. - rgba_map[:, :3])*(shad[:, num.newaxis])) 

170 elif blend_mode == 'multiply': 

171 rgba_map[:, :3] *= shad[:, num.newaxis] 

172 

173 return rgba_map 

174 

175 

176def plot_directivity( 

177 engine, source, store_id, 

178 distance=300*km, azi_begin=0., azi_end=360., dazi=1., 

179 phases={'P': 'first{stored:any_P}-10%', 

180 'S': 'last{stored:any_S}+50'}, 

181 interpolation='multilinear', 

182 target_depth=0.0, 

183 quantity='displacement', envelope=False, 

184 component='R', fmin=0.01, fmax=0.1, 

185 hillshade=True, cmap=None, 

186 plot_mt='full', show_phases=True, show_description=True, 

187 reverse_time=False, show_nucleations=True, axes=None, nthreads=0): 

188 ''' 

189 Plot the directivity and radiation characteristics of source models. 

190 

191 Synthetic seismic traces (R, T or Z) are forward-modelled at a defined 

192 radius, covering the full or partial azimuthal range and projected on a 

193 polar plot. Difference in the amplitude are enhanced by hillshading 

194 the data. 

195 

196 :param engine: Forward modelling engine 

197 :type engine: :py:class:`~pyrocko.gf.seismosizer.Engine` 

198 :param source: Parametrized source model 

199 :type source: :py:class:`~pyrocko.gf.seismosizer.Source` 

200 :param store_id: Store ID used for forward modelling 

201 :type store_id: str 

202 :param distance: Distance in [m] 

203 :type distance: float 

204 :param azi_begin: Begin azimuth in [deg] 

205 :type azi_begin: float 

206 :param azi_end: End azimuth in [deg] 

207 :type azi_end: float 

208 :param dazi: Delta azimuth, bin size [deg] 

209 :type dazi: float 

210 :param phases: Phases to define start and end of time window 

211 :type phases: :py:class:`dict` with :py:class:`str` keys and 

212 :py:class:`~pyrocko.gf.meta.Timing` values 

213 :param quantity: Seismogram quantity, default ``displacement`` 

214 :type quantity: str 

215 :param envelope: Plot envelope instead of seismic trace 

216 :type envelope: bool 

217 :param component: Forward modelled component, default ``R``. Choose from 

218 `RTZ` 

219 :type component: str 

220 :param fmin: Bandpass lower frequency [Hz], default ``0.01`` 

221 :type fmin: float 

222 :param fmax: Bandpass upper frequency [Hz], default ``0.1`` 

223 :type fmax: float 

224 :param hillshade: Enable hillshading, default ``True`` 

225 :type hillshade: bool 

226 :param cmap: Matplotlib colormap to use, default ``seismic``. 

227 When ``envelope`` is ``True`` the default colormap will be ``Reds``. 

228 :type cmap: str 

229 :param plot_mt: Plot a centered moment tensor, default ``full``. 

230 Choose from ``full, deviatoric, dc or False`` 

231 :type plot_mt: str, bool 

232 :param show_phases: Show annotations, default ``True`` 

233 :type show_phases: bool 

234 :param show_description: Show description, default ``True`` 

235 :type show_description: bool 

236 :param reverse_time: Reverse time axis. First phases arrive at the center, 

237 default ``False`` 

238 :type reverse_time: bool 

239 :param show_nucleations: Show nucleation piercing points on the moment 

240 tensor, default ``True`` 

241 :type show_nucleations: bool 

242 :param axes: Give axes to plot into 

243 :type axes: :py:class:`matplotlib.axes.Axes` 

244 :param nthreads: Number of threads used for forward modelling, 

245 default ``0`` - all available cores 

246 :type nthreads: int 

247 ''' 

248 

249 if axes is None: 

250 fig = plt.figure() 

251 ax = fig.add_subplot(111, polar=True) 

252 else: 

253 fig = axes.figure 

254 ax = axes 

255 

256 if envelope and cmap is None: 

257 cmap = 'Reds' 

258 elif cmap is None: 

259 cmap = 'seismic' 

260 

261 targets, azimuths = get_azimuthal_targets( 

262 store_id, source, distance, azi_begin, azi_end, dazi, 

263 depth=target_depth, 

264 interpolation=interpolation, 

265 components='R', quantity=quantity) 

266 ref_target = targets[0] 

267 store = engine.get_store(store_id) 

268 mt = source.pyrocko_moment_tensor(store=store, target=ref_target) 

269 

270 resp = engine.process(source, targets, nthreads=nthreads) 

271 data, times = get_seismogram_array( 

272 resp, fmin, fmax, 

273 component=component, envelope=envelope) 

274 

275 nucl_depth = source.depth 

276 nucl_distance = distance 

277 

278 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y'): 

279 try: 

280 iter(source.nucleation_x) 

281 nx = float(source.nucleation_x[0]) 

282 ny = float(source.nucleation_y[0]) 

283 

284 except TypeError: 

285 nx = source.nucleation_x 

286 ny = source.nucleation_y 

287 

288 nucl_distance += nx * source.length/2. 

289 nucl_depth += ny*num.sin(source.dip*d2r) * source.width/2. 

290 

291 if hasattr(source, 'anchor'): 

292 anch_x, anch_y = map_anchor[source.anchor] 

293 nucl_distance -= anch_x * source.length/2. 

294 nucl_depth -= anch_y*num.sin(source.dip*d2r) * source.width/2. 

295 

296 timings = [Timing(p) for p in phases.values()] 

297 phase_times = [store.t(t, source, ref_target) for t in timings] 

298 

299 tbegin = min(phase_times) 

300 tend = max(phase_times) 

301 tsel = num.logical_and(times >= tbegin, times <= tend) 

302 

303 data = data[:, tsel].T 

304 times = times[tsel] 

305 duration = times[-1] - times[0] 

306 

307 vmax = num.abs(data).max() 

308 cmw = ScalarMappable(cmap=cmap) 

309 cmw.set_array(data) 

310 cmw.set_clim(-vmax, vmax) 

311 

312 if envelope: 

313 cmw.set_clim(0., vmax) 

314 

315 ax.set_theta_zero_location('N') 

316 ax.set_theta_direction(-1) 

317 

318 strike_label = mt.strike1 

319 if hasattr(source, 'strike'): 

320 strike_label = source.strike 

321 

322 try: 

323 ax.set_rlabel_position(strike_label % 180. - 180.) 

324 except AttributeError: 

325 logger.warning('Old matplotlib version: cannot set label positions') 

326 

327 def r_fmt(v, p): 

328 if v < tbegin or v > tend: 

329 return '' 

330 return '%g s' % v 

331 

332 ax.yaxis.set_major_formatter(FuncFormatter(r_fmt)) 

333 if reverse_time: 

334 ax.set_rlim(times[0] - .3*duration, times[-1]) 

335 else: 

336 ax.set_rlim(times[-1] + .3*duration, times[0]) 

337 

338 ax.grid(zorder=20) 

339 

340 if isinstance(plot_mt, str): 

341 mt_size = .15 

342 beachball.plot_beachball_mpl( 

343 mt, ax, 

344 beachball_type=plot_mt, size=mt_size, 

345 size_units='axes', color_t=(0.7, 0.4, 0.4), 

346 position=(.5, .5), linewidth=1.) 

347 

348 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y')\ 

349 and show_nucleations: 

350 try: 

351 iter(source.nucleation_x) 

352 nucleation_x = source.nucleation_x 

353 nucleation_y = source.nucleation_y 

354 except TypeError: 

355 nucleation_x = [source.nucleation_x] 

356 nucleation_y = [source.nucleation_y] 

357 

358 for nx, ny in zip(nucleation_x, nucleation_y): 

359 angle = float(num.arctan2(ny, nx)) 

360 rtp = num.array([[1., angle, (90.-source.strike)*d2r]]) 

361 points = beachball.numpy_rtp2xyz(rtp) 

362 x, y = beachball.project(points, projection='lambert').T 

363 norm = num.sqrt(x**2 + y**2) 

364 x = x / norm * mt_size/2. 

365 y = y / norm * mt_size/2. 

366 ax.plot(x+.5, y+.5, 'x', ms=6, mew=2, mec='darkred', mfc='red', 

367 transform=ax.transAxes, zorder=10) 

368 

369 ax.grid(False) 

370 mesh = ax.pcolormesh( 

371 azimuths * d2r, times, data, 

372 cmap=cmw.cmap, norm=cmw.norm, shading='gouraud', zorder=0) 

373 

374 if hillshade: 

375 mesh.update_scalarmappable() 

376 color = mesh.get_facecolor() 

377 color = hillshade_seismogram_array( 

378 data, color, shad_lim=(.85, 1.), blend_mode='multiply') 

379 mesh.set_facecolor(color) 

380 

381 if show_phases: 

382 label_theta = 270. 

383 theta = num.linspace(0, 2*num.pi, 360) 

384 

385 for label, phase_str in phases.items(): 

386 phase = Timing(phase_str) 

387 

388 phase.offset = 0. 

389 phase.offset_is_slowness = False 

390 phase.offset_is_percent = False 

391 

392 time = store.t(phase, source, ref_target) 

393 times = num.full_like(theta, time) 

394 

395 ax.plot(theta, times, color='k', alpha=.3, lw=1., ls='--') 

396 

397 ax.text( 

398 label_theta*d2r, time, label, 

399 ha='left', color='k', fontsize='small') 

400 label_theta += 30. 

401 

402 if show_description: 

403 description = ( 

404 'Component {component:s}\n' 

405 'Distance {distance:g} km').format( 

406 component=component, distance=distance / km) 

407 

408 if fmin and fmax: 

409 description += '\nBandpass {fmin:g} - {fmax:g} Hz'.format( 

410 fmin=fmin, fmax=fmax) 

411 elif fmin: 

412 description += '\nHighpass {fmin:g} Hz'.format(fmin=fmin) 

413 elif fmax: 

414 description += '\nLowpass {fmax:g} Hz'.format(fmax=fmax) 

415 ax.text( 

416 -.05, -.05, description, 

417 fontsize='small', 

418 ha='left', va='bottom', transform=ax.transAxes) 

419 

420 cbar_label = QUANTITY_LABEL[quantity] 

421 if envelope: 

422 cbar_label = 'Envelope ' + cbar_label 

423 

424 cb = fig.colorbar( 

425 cmw, ax=ax, 

426 orientation='vertical', shrink=.8, pad=0.11) 

427 

428 cb.set_label(cbar_label) 

429 

430 if axes is None: 

431 plt.show() 

432 return resp 

433 

434 

435__all__ = ['plot_directivity'] 

436 

437 

438if __name__ == '__main__': 

439 engine = LocalEngine(store_superdirs=['.'], use_config=True) 

440 

441 rect_source = RectangularSource( 

442 depth=2.6*km, 

443 strike=240., 

444 dip=76.6, 

445 rake=-.4, 

446 anchor='top', 

447 

448 nucleation_x=-.57, 

449 nucleation_y=-.59, 

450 velocity=2070., 

451 

452 length=27*km, 

453 width=9.4*km, 

454 slip=1.4) 

455 

456 resp = plot_directivity( 

457 engine, rect_source, 'crust2_ib', 

458 dazi=5, component='R', quantity='displacement', envelope=True)