1# http://pyrocko.org - GPLv3
2#
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
103 if fmin is not None:
104 tr.highpass(4, fmin, demean=False)
106 if fmax is not None:
107 tr.lowpass(4, fmax, demean=False)
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)
113 for tr in traces:
114 tr.extend(tmin, tmax, fillmethod='repeat')
115 if envelope:
116 tr.envelope()
118 data = num.array([tr.get_ydata() for tr in traces])
119 nsamples = data.shape[1]
120 return data, num.linspace(tmin, tmax, nsamples)
123def hillshade(array, azimuth, angle_altitude):
124 azimuth = 360.0 - azimuth
125 azi = azimuth * r2d
126 alt = angle_altitude * r2d
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)
132 shaded = num.sin(alt)*num.sin(slope) \
133 + num.cos(alt)*num.cos(slope)*num.cos((azi - num.pi/2.) - aspect)
135 return (shaded + 1.)/2.
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
148 # convolution of two 2D arrays
149 shad = im_conv(seismogram_array, ramp.T).ravel()
150 shad *= -1.
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)
158 shad[shad > percentile98] = percentile98
159 shad[shad < percentile2] = percentile2
161 # # normalize shading
162 shad -= num.nanmin(shad)
163 shad /= num.nanmax(shad)
165 # # reduce range to balance gray color
166 shad *= shad_lim[1] - shad_lim[0]
167 shad += shad_lim[0]
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]
174 return rgba_map
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.
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.
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 '''
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
257 if envelope and cmap is None:
258 cmap = 'Reds'
259 elif cmap is None:
260 cmap = 'seismic'
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)
271 resp = engine.process(source, targets, nthreads=nthreads)
272 data, times = get_seismogram_array(
273 resp, fmin, fmax,
274 component=component, envelope=envelope)
276 nucl_depth = source.depth
277 nucl_distance = distance
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])
285 except TypeError:
286 nx = source.nucleation_x
287 ny = source.nucleation_y
289 nucl_distance += nx * source.length/2.
290 nucl_depth += ny*num.sin(source.dip*d2r) * source.width/2.
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.
297 timings = [Timing(p) for p in phases.values()]
298 phase_times = [store.t(t, source, ref_target) for t in timings]
300 tbegin = min(phase_times)
301 tend = max(phase_times)
302 tsel = num.logical_and(times >= tbegin, times <= tend)
304 data = data[:, tsel].T
305 times = times[tsel]
306 duration = times[-1] - times[0]
308 vmax = num.abs(data).max()
309 cmw = ScalarMappable(cmap=cmap)
310 cmw.set_array(data)
311 cmw.set_clim(-vmax, vmax)
313 if envelope:
314 cmw.set_clim(0., vmax)
316 ax.set_theta_zero_location('N')
317 ax.set_theta_direction(-1)
319 strike_label = mt.strike1
320 if hasattr(source, 'strike'):
321 strike_label = source.strike
323 try:
324 ax.set_rlabel_position(strike_label % 180. - 180.)
325 except AttributeError:
326 logger.warning('Old matplotlib version: cannot set label positions')
328 def r_fmt(v, p):
329 if v < tbegin or v > tend:
330 return ''
331 return '%g s' % v
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])
339 ax.grid(zorder=20)
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.)
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]
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)
370 mesh = ax.pcolormesh(
371 azimuths * d2r, times, data,
372 cmap=cmw.cmap, norm=cmw.norm, shading='gouraud', zorder=0)
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)
381 if show_phases:
382 label_theta = 270.
383 theta = num.linspace(0, 2*num.pi, 360)
385 for label, phase_str in phases.items():
386 phase = Timing(phase_str)
388 phase.offset = 0.
389 phase.offset_is_slowness = False
390 phase.offset_is_percent = False
392 time = store.t(phase, source, ref_target)
393 times = num_full_like(theta, time)
395 ax.plot(theta, times, color='k', alpha=.3, lw=1., ls='--')
397 ax.text(
398 label_theta*d2r, time, label,
399 ha='left', color='k', fontsize='small')
400 label_theta += 30.
402 if show_description:
403 description = (
404 'Component {component:s}\n'
405 'Distance {distance:g} km').format(
406 component=component, distance=distance / km)
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)
420 cbar_label = QUANTITY_LABEL[quantity]
421 if envelope:
422 cbar_label = 'Envelope ' + cbar_label
424 cb = fig.colorbar(
425 cmw, ax=ax,
426 orientation='vertical', shrink=.8, pad=0.11)
428 cb.set_label(cbar_label)
430 if axes is None:
431 plt.show()
432 return resp
435__all__ = ['plot_directivity']
438if __name__ == '__main__':
439 engine = LocalEngine(store_superdirs=['.'], use_config=True)
441 rect_source = RectangularSource(
442 depth=2.6*km,
443 strike=240.,
444 dip=76.6,
445 rake=-.4,
446 anchor='top',
448 nucleation_x=-.57,
449 nucleation_y=-.59,
450 velocity=2070.,
452 length=27*km,
453 width=9.4*km,
454 slip=1.4)
456 resp = plot_directivity(
457 engine, rect_source, 'crust2_ib',
458 dazi=5, component='R', quantity='displacement', envelope=True)