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 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y'):
277 try:
278 iter(source.nucleation_x)
279 nx = float(source.nucleation_x[0])
280 ny = float(source.nucleation_y[0])
282 except TypeError:
283 nx = source.nucleation_x
284 ny = source.nucleation_y
286 nucl_distance += nx * source.length/2.
287 nucl_depth += ny*num.sin(source.dip*d2r) * source.width/2.
289 if hasattr(source, 'anchor'):
290 anch_x, anch_y = map_anchor[source.anchor]
291 nucl_distance -= anch_x * source.length/2.
292 nucl_depth -= anch_y*num.sin(source.dip*d2r) * source.width/2.
294 timings = [Timing(p) for p in phases.values()]
295 phase_times = [store.t(t, source, ref_target) for t in timings]
297 tbegin = min(phase_times)
298 tend = max(phase_times)
299 tsel = num.logical_and(times >= tbegin, times <= tend)
301 data = data[:, tsel].T
302 times = times[tsel]
303 duration = times[-1] - times[0]
305 vmax = num.abs(data).max()
306 cmw = ScalarMappable(cmap=cmap)
307 cmw.set_array(data)
308 cmw.set_clim(-vmax, vmax)
310 if envelope:
311 cmw.set_clim(0., vmax)
313 ax.set_theta_zero_location('N')
314 ax.set_theta_direction(-1)
316 strike_label = mt.strike1
317 if hasattr(source, 'strike'):
318 strike_label = source.strike
320 try:
321 ax.set_rlabel_position(strike_label % 180.)
322 except AttributeError:
323 logger.warn('Old matplotlib version: cannot set label positions')
325 def r_fmt(v, p):
326 if v < tbegin or v > tend:
327 return ''
328 return '%g s' % v
330 ax.yaxis.set_major_formatter(FuncFormatter(r_fmt))
331 if reverse_time:
332 ax.set_rlim(times[0] - .3*duration, times[-1])
333 else:
334 ax.set_rlim(times[-1] + .3*duration, times[0])
336 ax.grid(zorder=20)
338 if isinstance(plot_mt, str):
339 mt_size = .15
340 beachball.plot_beachball_mpl(
341 mt, ax,
342 beachball_type=plot_mt, size=mt_size,
343 size_units='axes', color_t=(0.7, 0.4, 0.4),
344 position=(.5, .5), linewidth=1.)
346 if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y')\
347 and show_nucleations:
348 try:
349 iter(source.nucleation_x)
350 nucleation_x = source.nucleation_x
351 nucleation_y = source.nucleation_y
352 except TypeError:
353 nucleation_x = [source.nucleation_x]
354 nucleation_y = [source.nucleation_y]
356 for nx, ny in zip(nucleation_x, nucleation_y):
357 angle = float(num.arctan2(ny, nx))
358 rtp = num.array([[1., angle, (90.-source.strike)*d2r]])
359 points = beachball.numpy_rtp2xyz(rtp)
360 x, y = beachball.project(points, projection='lambert').T
361 norm = num.sqrt(x**2 + y**2)
362 x = x / norm * mt_size/2.
363 y = y / norm * mt_size/2.
364 ax.plot(x+.5, y+.5, 'x', ms=6, mew=2, mec='darkred', mfc='red',
365 transform=ax.transAxes, zorder=10)
367 mesh = ax.pcolormesh(
368 azimuths * d2r, times, data,
369 cmap=cmw.cmap, shading='gouraud', zorder=0)
371 if hillshade:
372 mesh.update_scalarmappable()
373 color = mesh.get_facecolor()
374 color = hillshade_seismogram_array(
375 data, color, shad_lim=(.85, 1.), blend_mode='multiply')
376 mesh.set_facecolor(color)
378 if show_phases:
379 label_theta = 270.
380 theta = num.linspace(0, 2*num.pi, 360)
382 for label, phase_str in phases.items():
383 phase = Timing(phase_str)
385 phase.offset = 0.
386 phase.offset_is_slowness = False
387 phase.offset_is_percent = False
389 time = store.t(phase, source, ref_target)
390 times = num_full_like(theta, time)
392 ax.plot(theta, times, color='k', alpha=.3, lw=1., ls='--')
394 ax.text(
395 label_theta*d2r, time, label,
396 ha='left', color='k', fontsize='small')
397 label_theta += 30.
399 if show_description:
400 description = (
401 'Component {component:s}\n'
402 'Distance {distance:g} km').format(
403 component=component, distance=distance / km)
405 if fmin and fmax:
406 description += '\nBandpass {fmin:g} - {fmax:g} Hz'.format(
407 fmin=fmin, fmax=fmax)
408 elif fmin:
409 description += '\nHighpass {fmin:g} Hz'.format(fmin=fmin)
410 elif fmax:
411 description += '\nLowpass {fmax:g} Hz'.format(fmax=fmax)
412 ax.text(
413 -.05, -.05, description,
414 fontsize='small',
415 ha='left', va='bottom', transform=ax.transAxes)
417 cbar_label = QUANTITY_LABEL[quantity]
418 if envelope:
419 cbar_label = 'Envelope ' + cbar_label
421 cb = fig.colorbar(
422 cmw, ax=ax,
423 orientation='vertical', shrink=.8, pad=0.11)
425 cb.set_label(cbar_label)
427 if axes is None:
428 plt.show()
429 return resp
432__all__ = ['plot_directivity']
435if __name__ == '__main__':
436 engine = LocalEngine(store_superdirs=['.'], use_config=True)
438 rect_source = RectangularSource(
439 depth=2.6*km,
440 strike=240.,
441 dip=76.6,
442 rake=-.4,
443 anchor='top',
445 nucleation_x=-.57,
446 nucleation_y=-.59,
447 velocity=2070.,
449 length=27*km,
450 width=9.4*km,
451 slip=1.4)
453 resp = plot_directivity(
454 engine, rect_source, 'crust2_ib',
455 dazi=5, component='R', quantity='displacement', envelope=True)