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

16from pyrocko.util import num_full_like

19km = 1e3

20r2d = 180. / num.pi

21d2r = num.pi / 180.

23logger = logging.getLogger(__name__)

26QUANTITY_LABEL = {

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

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

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

30}

33def get_azimuthal_targets(

34 store_id, source, radius,

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

36 depth=0.0,

37 interpolation='multilinear',

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

40 assert dazi > 0.

41 assert azi_begin < azi_end

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

44 assert nstations > 0

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

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

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

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

51 coords *= radius

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

54 for comp in components:

55 assert comp in dips.keys()

57 target_kwargs = dict(

58 quantity=quantity,

59 interpolation=interpolation,

60 store_id=store_id)

62 targets = [

63 Target(

64 lat=source.lat,

65 lon=source.lon,

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

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

68 depth=depth,

69 azimuth={

70 'R': azi,

71 'T': azi+90.,

72 'Z': 0.

73 }[channel],

74 dip=dips[channel],

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

76 **target_kwargs)

77 for iazi, azi in enumerate(azimuths)

78 for channel in components]

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

81 target.azimuth = azi

82 target.dazi = dazi

84 return targets, azimuths

87def get_seismogram_array(

88 response, fmin=None, fmax=None,

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

90 resp = response

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

93 tmin = None

94 tmax = None

95 traces = []

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

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

99 continue

100 assert hasattr(target, 'azimuth')

101 assert target.dazi

102

103 if fmin is not None:

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

105

106 if fmax is not None:

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

108

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

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

111 traces.append(tr)

112

113 for tr in traces:

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

115 if envelope:

116 tr.envelope()

117

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

119 nsamples = data.shape[1]

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

121

122

123def hillshade(array, azimuth, angle_altitude):

124 azimuth = 360.0 - azimuth

125 azi = azimuth * r2d

126 alt = angle_altitude * r2d

127

128 x, y = num.gradient(array)

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

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

131

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

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

134

135 return (shaded + 1.)/2.

136

137

138def hillshade_seismogram_array(

139 seismogram_array, rgba_map,

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

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

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

143 from scipy.ndimage import convolve as im_conv

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

145 # from upper left

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

147

148 # convolution of two 2D arrays

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

150 shad *= -1.

151

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

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

154 # # shades helps

155 percentile2 = num.percentile(shad, 2.0)

156 percentile98 = num.percentile(shad, 98.0)

157

158 shad[shad > percentile98] = percentile98

159 shad[shad < percentile2] = percentile2

160

161 # # normalize shading

162 shad -= num.nanmin(shad)

163 shad /= num.nanmax(shad)

164

165 # # reduce range to balance gray color

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

167 shad += shad_lim[0]

168

169 if blend_mode == 'screen':

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

171 elif blend_mode == 'multiply':

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

173

174 return rgba_map

175

176

177def plot_directivity(

178 engine, source, store_id,

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

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

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

182 interpolation='multilinear',

183 target_depth=0.0,

184 quantity='displacement', envelope=False,

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

186 hillshade=True, cmap=None,

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

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

189 '''

190 Plot the directivity and radiation characteristics of source models.

191

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

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

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

195 the data.

196

197 :param engine: Forward modelling engine

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

199 :param source: Parametrized source model

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

201 :param store_id: Store ID used for forward modelling

202 :type store_id: str

203 :param distance: Distance in [m]

204 :type distance: float

205 :param azi_begin: Begin azimuth in [deg]

206 :type azi_begin: float

207 :param azi_end: End azimuth in [deg]

208 :type azi_end: float

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

210 :type dazi: float

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

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

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

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

215 :type quantity: str

216 :param envelope: Plot envelope instead of seismic trace

217 :type envelope: bool

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

219 `RTZ`

220 :type component: str

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

222 :type fmin: float

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

224 :type fmax: float

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

226 :type hillshade: bool

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

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

229 :type cmap: str

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

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

232 :type plot_mt: str, bool

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

234 :type show_phases: bool

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

236 :type show_description: bool

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

238 default ``False``

239 :type reverse_time: bool

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

241 tensor, default ``True``

242 :type show_nucleations: bool

243 :param axes: Give axes to plot into

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

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

246 default ``0`` - all available cores

247 :type nthreads: int

248 '''

249

250 if axes is None:

251 fig = plt.figure()

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

253 else:

254 fig = axes.figure

255 ax = axes

256

257 if envelope and cmap is None:

258 cmap = 'Reds'

259 elif cmap is None:

260 cmap = 'seismic'

261

262 targets, azimuths = get_azimuthal_targets(

263 store_id, source, distance, azi_begin, azi_end, dazi,

264 depth=target_depth,

265 interpolation=interpolation,

266 components='R', quantity=quantity)

267 ref_target = targets[0]

268 store = engine.get_store(store_id)

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

270

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

272 data, times = get_seismogram_array(

273 resp, fmin, fmax,

274 component=component, envelope=envelope)

275

276 nucl_depth = source.depth

277 nucl_distance = distance

278

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

280 try:

281 iter(source.nucleation_x)

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

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

284

285 except TypeError:

286 nx = source.nucleation_x

287 ny = source.nucleation_y

288

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

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

291

292 if hasattr(source, 'anchor'):

293 anch_x, anch_y = map_anchor[source.anchor]

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

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

296

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

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

299

300 tbegin = min(phase_times)

301 tend = max(phase_times)

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

303

304 data = data[:, tsel].T

305 times = times[tsel]

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

307

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

309 cmw = ScalarMappable(cmap=cmap)

310 cmw.set_array(data)

311 cmw.set_clim(-vmax, vmax)

312

313 if envelope:

314 cmw.set_clim(0., vmax)

315

316 ax.set_theta_zero_location('N')

317 ax.set_theta_direction(-1)

318

319 strike_label = mt.strike1

320 if hasattr(source, 'strike'):

321 strike_label = source.strike

322

323 try:

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

325 except AttributeError:

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

327

328 def r_fmt(v, p):

329 if v < tbegin or v > tend:

330 return ''

331 return '%g s' % v

332

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

334 if reverse_time:

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

336 else:

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

338

339 ax.grid(zorder=20)

340

341 if isinstance(plot_mt, str):

342 mt_size = .15

343 beachball.plot_beachball_mpl(

344 mt, ax,

345 beachball_type=plot_mt, size=mt_size,

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

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

348

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

350 and show_nucleations:

351 try:

352 iter(source.nucleation_x)

353 nucleation_x = source.nucleation_x

354 nucleation_y = source.nucleation_y

355 except TypeError:

356 nucleation_x = [source.nucleation_x]

357 nucleation_y = [source.nucleation_y]

358

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

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

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

362 points = beachball.numpy_rtp2xyz(rtp)

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

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

365 x = x / norm * mt_size/2.

366 y = y / norm * mt_size/2.

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

368 transform=ax.transAxes, zorder=10)

369

370 ax.grid(False)

371 mesh = ax.pcolormesh(

372 azimuths * d2r, times, data,

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

374

375 if hillshade:

376 mesh.update_scalarmappable()

377 color = mesh.get_facecolor()

378 color = hillshade_seismogram_array(

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

380 mesh.set_facecolor(color)

381

382 if show_phases:

383 label_theta = 270.

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

385

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

387 phase = Timing(phase_str)

388

389 phase.offset = 0.

390 phase.offset_is_slowness = False

391 phase.offset_is_percent = False

392

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

394 times = num_full_like(theta, time)

395

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

397

398 ax.text(

399 label_theta*d2r, time, label,

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

401 label_theta += 30.

402

403 if show_description:

404 description = (

405 'Component {component:s}\n'

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

407 component=component, distance=distance / km)

408

409 if fmin and fmax:

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

411 fmin=fmin, fmax=fmax)

412 elif fmin:

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

414 elif fmax:

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

416 ax.text(

417 -.05, -.05, description,

418 fontsize='small',

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

420

421 cbar_label = QUANTITY_LABEL[quantity]

422 if envelope:

423 cbar_label = 'Envelope ' + cbar_label

424

425 cb = fig.colorbar(

426 cmw, ax=ax,

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

428

429 cb.set_label(cbar_label)

430

431 if axes is None:

432 plt.show()

433 return resp

434

435

436__all__ = ['plot_directivity']

437

438

439if __name__ == '__main__':

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

441

442 rect_source = RectangularSource(

443 depth=2.6*km,

444 strike=240.,

445 dip=76.6,

446 rake=-.4,

447 anchor='top',

448

449 nucleation_x=-.57,

450 nucleation_y=-.59,

451 velocity=2070.,

452

453 length=27*km,

454 width=9.4*km,

455 slip=1.4)

456

457 resp = plot_directivity(

458 engine, rect_source, 'crust2_ib',

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