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