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.abshilbert()
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 phase_begin: Start time of the window
212 :type phase_begin: :py:class:`~pyrocko.gf.meta.Timing`
213 :param phase_end: End time of the window
214 :type phase_end: :py:class:`~pyrocko.gf.meta.Timing`
215 :param target_depth: Receiver depth [m]
216 :type target_depth: float
217 :param quantity: Seismogram quantity, default ``displacement``
218 :type quantity: str
219 :param envelope: Plot envelope instead of seismic trace
220 :type envelope: bool
221 :param component: Forward modelled component, default ``R``. Choose from
222 `RTZ`
223 :type component: str
224 :param fmin: Bandpass lower frequency [Hz], default ``0.01``
225 :type fmin: float
226 :param fmax: Bandpass upper frequency [Hz], default ``0.1``
227 :type fmax: float
228 :param hillshade: Enable hillshading, default ``True``
229 :type hillshade: bool
230 :param cmap: Matplotlib colormap to use, default ``seismic``.
231 When ``envelope`` is ``True`` the default colormap will be ``Reds``.
232 :type cmap: str
233 :param plot_mt: Plot a centered moment tensor, default ``full``.
234 Choose from ``full, deviatoric, dc or False``
235 :type plot_mt: str, bool
236 :param show_phases: Show annotations, default ``True``
237 :type show_phases: bool
238 :param show_description: Show description, default ``True``
239 :type show_description: bool
240 :param reverse_time: Reverse time axis. First phases arrive at the center,
241 default ``False``
242 :type reverse_time: bool
243 :param show_nucleations: Show nucleation piercing points on the moment
244 tensor, default ``True``
245 :type show_nucleations: bool
246 :param axes: Give axes to plot into
247 :type axes: :py:class:`matplotlib.axes.Axes`
248 :param nthreads: Number of threads used for forward modelling,
249 default ``0`` - all available cores
250 :type nthreads: int
251 '''
253 if axes is None:
254 fig = plt.figure()
255 ax = fig.add_subplot(111, polar=True)
256 else:
257 fig = axes.figure
258 ax = axes
260 if envelope and cmap is None:
261 cmap = 'Reds'
262 elif cmap is None:
263 cmap = 'seismic'
265 targets, azimuths = get_azimuthal_targets(
266 store_id, source, distance, azi_begin, azi_end, dazi,
267 depth=target_depth,
268 interpolation=interpolation,
269 components='R', quantity=quantity)
270 ref_target = targets[0]
271 store = engine.get_store(store_id)
272 mt = source.pyrocko_moment_tensor(store=store, target=ref_target)
274 resp = engine.process(source, targets, nthreads=nthreads)
275 data, times = get_seismogram_array(
276 resp, fmin, fmax,
277 component=component, envelope=envelope)
279 nucl_depth = source.depth
280 nucl_distance = distance
282 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y'):
283 try:
284 iter(source.nucleation_x)
285 nx = float(source.nucleation_x[0])
286 ny = float(source.nucleation_y[0])
288 except TypeError:
289 nx = source.nucleation_x
290 ny = source.nucleation_y
292 nucl_distance += nx * source.length/2.
293 nucl_depth += ny*num.sin(source.dip*d2r) * source.width/2.
295 if hasattr(source, 'anchor'):
296 anch_x, anch_y = map_anchor[source.anchor]
297 nucl_distance -= anch_x * source.length/2.
298 nucl_depth -= anch_y*num.sin(source.dip*d2r) * source.width/2.
300 timings = [Timing(p) for p in phases.values()]
301 phase_times = [store.t(t, source, ref_target) for t in timings]
303 tbegin = min(phase_times)
304 tend = max(phase_times)
305 tsel = num.logical_and(times >= tbegin, times <= tend)
307 data = data[:, tsel].T
308 times = times[tsel]
309 duration = times[-1] - times[0]
311 vmax = num.abs(data).max()
312 cmw = ScalarMappable(cmap=cmap)
313 cmw.set_array(data)
314 cmw.set_clim(-vmax, vmax)
316 if envelope:
317 cmw.set_clim(0., vmax)
319 ax.set_theta_zero_location('N')
320 ax.set_theta_direction(-1)
322 strike_label = mt.strike1
323 if hasattr(source, 'strike'):
324 strike_label = source.strike
326 try:
327 ax.set_rlabel_position(strike_label % 180.)
328 except AttributeError:
329 logger.warn('Old matplotlib version: cannot set label positions')
331 def r_fmt(v, p):
332 if v < tbegin or v > tend:
333 return ''
334 return '%g s' % v
336 ax.yaxis.set_major_formatter(FuncFormatter(r_fmt))
337 if reverse_time:
338 ax.set_rlim(times[0] - .3*duration, times[-1])
339 else:
340 ax.set_rlim(times[-1] + .3*duration, times[0])
342 ax.grid(zorder=20)
344 if isinstance(plot_mt, str):
345 mt_size = .15
346 beachball.plot_beachball_mpl(
347 mt, ax,
348 beachball_type=plot_mt, size=mt_size,
349 size_units='axes', color_t=(0.7, 0.4, 0.4),
350 position=(.5, .5), linewidth=1.)
352 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y')\
353 and show_nucleations:
354 try:
355 iter(source.nucleation_x)
356 nucleation_x = source.nucleation_x
357 nucleation_y = source.nucleation_y
358 except TypeError:
359 nucleation_x = [source.nucleation_x]
360 nucleation_y = [source.nucleation_y]
362 for nx, ny in zip(nucleation_x, nucleation_y):
363 angle = float(num.arctan2(ny, nx))
364 rtp = num.array([[1., angle, (90.-source.strike)*d2r]])
365 points = beachball.numpy_rtp2xyz(rtp)
366 x, y = beachball.project(points, projection='lambert').T
367 norm = num.sqrt(x**2 + y**2)
368 x = x / norm * mt_size/2.
369 y = y / norm * mt_size/2.
370 ax.plot(x+.5, y+.5, 'x', ms=6, mew=2, mec='darkred', mfc='red',
371 transform=ax.transAxes, zorder=10)
373 mesh = ax.pcolormesh(
374 azimuths * d2r, times, data,
375 cmap=cmw.cmap, norm=cmw.norm, shading='gouraud', zorder=0)
377 if hillshade:
378 mesh.update_scalarmappable()
379 color = mesh.get_facecolor()
380 color = hillshade_seismogram_array(
381 data, color, shad_lim=(.85, 1.), blend_mode='multiply')
382 mesh.set_facecolor(color)
384 if show_phases:
385 label_theta = 270.
386 theta = num.linspace(0, 2*num.pi, 360)
388 for label, phase_str in phases.items():
389 phase = Timing(phase_str)
391 phase.offset = 0.
392 phase.offset_is_slowness = False
393 phase.offset_is_percent = False
395 time = store.t(phase, source, ref_target)
396 times = num_full_like(theta, time)
398 ax.plot(theta, times, color='k', alpha=.3, lw=1., ls='--')
400 ax.text(
401 label_theta*d2r, time, label,
402 ha='left', color='k', fontsize='small')
403 label_theta += 30.
405 if show_description:
406 description = (
407 'Component {component:s}\n'
408 'Distance {distance:g} km').format(
409 component=component, distance=distance / km)
411 if fmin and fmax:
412 description += '\nBandpass {fmin:g} - {fmax:g} Hz'.format(
413 fmin=fmin, fmax=fmax)
414 elif fmin:
415 description += '\nHighpass {fmin:g} Hz'.format(fmin=fmin)
416 elif fmax:
417 description += '\nLowpass {fmax:g} Hz'.format(fmax=fmax)
418 ax.text(
419 -.05, -.05, description,
420 fontsize='small',
421 ha='left', va='bottom', transform=ax.transAxes)
423 cbar_label = QUANTITY_LABEL[quantity]
424 if envelope:
425 cbar_label = 'Envelope ' + cbar_label
427 cb = fig.colorbar(
428 cmw, ax=ax,
429 orientation='vertical', shrink=.8, pad=0.11)
431 cb.set_label(cbar_label)
433 if axes is None:
434 plt.show()
435 return resp
438__all__ = ['plot_directivity']
441if __name__ == '__main__':
442 engine = LocalEngine(store_superdirs=['.'], use_config=True)
444 rect_source = RectangularSource(
445 depth=2.6*km,
446 strike=240.,
447 dip=76.6,
448 rake=-.4,
449 anchor='top',
451 nucleation_x=-.57,
452 nucleation_y=-.59,
453 velocity=2070.,
455 length=27*km,
456 width=9.4*km,
457 slip=1.4)
459 resp = plot_directivity(
460 engine, rect_source, 'crust2_ib',
461 dazi=5, component='R', quantity='displacement', envelope=True)