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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

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

6import logging

7import numpy as num

8import matplotlib.pyplot as plt

10from matplotlib.cm import ScalarMappable

11from matplotlib.ticker import FuncFormatter

13from pyrocko.plot import beachball

14from pyrocko.gf.meta import Timing

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

18km = 1e3

19r2d = 180. / num.pi

20d2r = num.pi / 180.

22logger = logging.getLogger(__name__)

25QUANTITY_LABEL = {

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

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

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

29}

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'):

39 assert dazi > 0.

40 assert azi_begin < azi_end

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

43 assert nstations > 0

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

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

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

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

50 coords *= radius

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

53 for comp in components:

54 assert comp in dips.keys()

56 target_kwargs = dict(

57 quantity=quantity,

58 interpolation=interpolation,

59 store_id=store_id)

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]

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

80 target.azimuth = azi

81 target.dazi = dazi

83 return targets, azimuths

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'

92 tmin = None

93 tmax = None

94 traces = []

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)