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 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='displacement',
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 azimuth={
68 'R': azi,
69 'T': azi+90.,
70 'Z': 0.
71 }[channel],
72 dip=dips[channel],
73 codes=('', 'S%01d' % iazi, '', channel),
74 **target_kwargs)
75 for iazi, azi in enumerate(azimuths)
76 for channel in components]
78 for target, azi in zip(targets, azimuths):
79 target.azimuth = azi
80 target.dazi = dazi
82 return targets, azimuths
85def get_seismogram_array(
86 response, fmin=None, fmax=None,
87 component='R', envelope=False):
88 resp = response
89 assert len(resp.request.sources) == 1, 'more than one source in response'
91 tmin = None
92 tmax = None
93 traces = []
95 for _, target, tr in response.iter_results():
96 if target.codes[-1] != component:
97 continue
98 assert hasattr(target, 'azimuth')
99 assert target.dazi
101 if fmin and fmax:
102 tr.bandpass(2, fmin, fmax)
103 elif fmin:
104 tr.highpass(4, fmin)
105 elif fmax:
106 tr.lowpass(4, fmax)
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)
112 for tr in traces:
113 tr.extend(tmin, tmax, fillmethod='repeat')
114 if envelope:
115 tr.abshilbert()
117 data = num.array([tr.get_ydata() for tr in traces])
118 data -= data.mean()
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 quantity='displacement', envelope=False,
183 component='R', fmin=0.01, fmax=0.1,
184 hillshade=True, cmap=None,
185 plot_mt='full', show_phases=True, show_description=True,
186 reverse_time=False, show_nucleations=True, axes=None, nthreads=0):
187 '''
188 Plot the directivity and radiation characteristics of source models.
190 Synthetic seismic traces (R, T or Z) are forward-modelled at a defined
191 radius, covering the full or partial azimuthal range and projected on a
192 polar plot. Difference in the amplitude are enhanced by hillshading
193 the data.
195 :param engine: Forward modelling engine
196 :type engine: :py:class:`~pyrocko.gf.seismosizer.Engine`
197 :param source: Parametrized source model
198 :type source: :py:class:`~pyrocko.gf.seismosizer.Source`
199 :param store_id: Store ID used for forward modelling
200 :type store_id: str
201 :param distance: Distance in [m]
202 :type distance: float
203 :param azi_begin: Begin azimuth in [deg]
204 :type azi_begin: float
205 :param azi_end: End azimuth in [deg]
206 :type azi_end: float
207 :param dazi: Delta azimuth, bin size [deg]
208 :type dazi: float
209 :param phase_begin: Start time of the window
210 :type phase_begin: :py:class:`~pyrocko.gf.meta.Timing`
211 :param phase_end: End time of the window
212 :type phase_end: :py:class:`~pyrocko.gf.meta.Timing`
213 :param quantity: Seismogram quantity, default ``displacement``
214 :type quantity: str
215 :param envelope: Plot envelop 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: Matplotlit 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 desciption, 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 '''
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
256 if envelope and cmap is None:
257 cmap = 'Reds'
258 elif cmap is None:
259 cmap = 'seismic'
261 targets, azimuths = get_azimuthal_targets(
262 store_id, source, distance, azi_begin, azi_end, dazi,
263 components='R', quantity=quantity)
264 ref_target = targets[0]
265 store = engine.get_store(store_id)
266 mt = source.pyrocko_moment_tensor(store=store, target=ref_target)
268 resp = engine.process(source, targets, nthreads=nthreads)
269 data, times = get_seismogram_array(
270 resp, fmin, fmax,
271 component=component, envelope=envelope)
273 nucl_depth = source.depth
274 nucl_distance = distance
276 anch_x, anch_y = map_anchor[source.anchor]
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])
284 except TypeError:
285 nx = source.nucleation_x
286 ny = source.nucleation_y
288 nucl_distance += nx * source.length/2.
289 nucl_depth += ny*num.sin(source.dip*d2r) * source.width/2.
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.
296 timings = [Timing(p) for p in phases.values()]
297 phase_times = [store.t(t, source, ref_target) for t in timings]
299 tbegin = min(phase_times)
300 tend = max(phase_times)
301 tsel = num.logical_and(times >= tbegin, times <= tend)
303 data = data[:, tsel].T
304 times = times[tsel]
305 duration = times[-1] - times[0]
307 vmax = num.abs(data).max()
308 cmw = ScalarMappable(cmap=cmap)
309 cmw.set_array(data)
310 cmw.set_clim(-vmax, vmax)
312 if envelope:
313 cmw.set_clim(0., vmax)
315 ax.set_theta_zero_location('N')
316 ax.set_theta_direction(-1)
318 strike_label = mt.strike1
319 if hasattr(source, 'strike'):
320 strike_label = source.strike
322 try:
323 ax.set_rlabel_position(strike_label % 180. - 180.)
324 except AttributeError:
325 logger.warn('Old matplotlib version: cannot set label positions')
327 def r_fmt(v, p):
328 if v < tbegin or v > tend:
329 return ''
330 return '%g s' % v
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])
338 ax.grid(zorder=20)
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.)
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]
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)
369 mesh = ax.pcolormesh(
370 azimuths * d2r, times, data,
371 cmap=cmw.cmap, shading='gouraud', zorder=0)
373 if hillshade:
374 mesh.update_scalarmappable()
375 color = mesh.get_facecolor()
376 color = hillshade_seismogram_array(
377 data, color, shad_lim=(.85, 1.), blend_mode='multiply')
378 mesh.set_facecolor(color)
380 if show_phases:
381 label_theta = 270.
382 theta = num.linspace(0, 2*num.pi, 360)
384 for label, phase_str in phases.items():
385 phase = Timing(phase_str)
387 phase.offset = 0.
388 phase.offset_is_slowness = False
389 phase.offset_is_percent = False
391 time = store.t(phase, source, ref_target)
392 times = num_full_like(theta, time)
394 ax.plot(theta, times, color='k', alpha=.3, lw=1., ls='--')
396 ax.text(
397 label_theta*d2r, time, label,
398 ha='left', color='k', fontsize='small')
399 label_theta += 30.
401 if show_description:
402 description = (
403 'Component {component:s}\n'
404 'Distance {distance:g} km').format(
405 component=component, distance=distance / km)
407 if fmin and fmax:
408 description += '\nBandpass {fmin:g} - {fmax:g} Hz'.format(
409 fmin=fmin, fmax=fmax)
410 elif fmin:
411 description += '\nHighpass {fmin:g} Hz'.format(fmin=fmin)
412 elif fmax:
413 description += '\nLowpass {fmax:g} Hz'.format(fmax=fmax)
414 ax.text(
415 -.05, -.05, description,
416 fontsize='small',
417 ha='left', va='bottom', transform=ax.transAxes)
419 cbar_label = QUANTITY_LABEL[quantity]
420 if envelope:
421 cbar_label = 'Envelope ' + cbar_label
423 cb = fig.colorbar(
424 cmw, ax=ax,
425 orientation='vertical', shrink=.8, pad=0.11)
427 cb.set_label(cbar_label)
429 if axes is None:
430 plt.show()
431 return resp
434__all__ = ['plot_directivity']
437if __name__ == '__main__':
438 engine = LocalEngine(store_superdirs=['.'], use_config=True)
440 rect_source = RectangularSource(
441 depth=2.6*km,
442 strike=240.,
443 dip=76.6,
444 rake=-.4,
445 anchor='top',
447 nucleation_x=-.57,
448 nucleation_y=-.59,
449 velocity=2070.,
451 length=27*km,
452 width=9.4*km,
453 slip=1.4)
455 resp = plot_directivity(
456 engine, rect_source, 'crust2_ib',
457 dazi=5, component='R', quantity='displacement', envelope=True)