Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/directivity.py: 87%
208 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-01-15 12:05 +0000
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
18km = 1e3
19r2d = 180. / num.pi
20d2r = num.pi / 180.
22logger = logging.getLogger(__name__)
25QUANTITY_LABEL = {
26 'displacement': 'Displacement [m]',
27 'velocity': 'Velocity [m/s]',
28 'acceleration': 'Acceleration [m/s²]'
29}
32def get_azimuthal_targets(
33 store_id, source, radius,
34 azi_begin=0., azi_end=360., dazi=1.,
35 depth=0.0,
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=quantity,
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 depth=depth,
68 azimuth={
69 'R': azi,
70 'T': azi+90.,
71 'Z': 0.
72 }[channel],
73 dip=dips[channel],
74 codes=('', 'S%01d' % iazi, '', channel),
75 **target_kwargs)
76 for iazi, azi in enumerate(azimuths)
77 for channel in components]
79 for target, azi in zip(targets, azimuths):
80 target.azimuth = azi
81 target.dazi = dazi
83 return targets, azimuths
86def get_seismogram_array(
87 response, fmin=None, fmax=None,
88 component='R', envelope=False):
89 resp = response
90 assert len(resp.request.sources) == 1, 'more than one source in response'
92 tmin = None
93 tmax = None
94 traces = []
96 for _, target, tr in response.iter_results():
97 if target.codes[-1] != component:
98 continue
99 assert hasattr(target, 'azimuth')
100 assert target.dazi
102 if fmin is not None:
103 tr.highpass(4, fmin, demean=False)
105 if fmax is not None:
106 tr.lowpass(4, fmax, demean=False)
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.envelope()
117 data = num.array([tr.get_ydata() for tr in traces])
118 nsamples = data.shape[1]
119 return data, num.linspace(tmin, tmax, nsamples)
122def hillshade(array, azimuth, angle_altitude):
123 azimuth = 360.0 - azimuth
124 azi = azimuth * r2d
125 alt = angle_altitude * r2d
127 x, y = num.gradient(array)
128 slope = num.pi/2. - num.arctan(num.sqrt(x*x + y*y))
129 aspect = num.arctan2(-x, y)
131 shaded = num.sin(alt)*num.sin(slope) \
132 + num.cos(alt)*num.cos(slope)*num.cos((azi - num.pi/2.) - aspect)
134 return (shaded + 1.)/2.
137def hillshade_seismogram_array(
138 seismogram_array, rgba_map,
139 shad_lim=(.4, .98), contrast=1., blend_mode='multiply'):
140 assert blend_mode in ('multiply', 'screen'), 'unknown blend mode'
141 assert shad_lim[0] < shad_lim[1], 'bad shading limits'
142 from scipy.ndimage import convolve as im_conv
143 # Light source from somewhere above - psychologically the best choice
144 # from upper left
145 ramp = num.array([[1., 0.], [0., -1.]]) * contrast
147 # convolution of two 2D arrays
148 shad = im_conv(seismogram_array, ramp.T).ravel()
149 shad *= -1.
151 # if there are strong artifical edges in the data, shades get
152 # dominated by them. Cutting off the largest and smallest 2% of
153 # # shades helps
154 percentile2 = num.percentile(shad, 2.0)
155 percentile98 = num.percentile(shad, 98.0)
157 shad[shad > percentile98] = percentile98
158 shad[shad < percentile2] = percentile2
160 # # normalize shading
161 shad -= num.nanmin(shad)
162 shad /= num.nanmax(shad)
164 # # reduce range to balance gray color
165 shad *= shad_lim[1] - shad_lim[0]
166 shad += shad_lim[0]
168 if blend_mode == 'screen':
169 rgba_map[:, :3] = 1. - ((1. - rgba_map[:, :3])*(shad[:, num.newaxis]))
170 elif blend_mode == 'multiply':
171 rgba_map[:, :3] *= shad[:, num.newaxis]
173 return rgba_map
176def plot_directivity(
177 engine, source, store_id,
178 distance=300*km, azi_begin=0., azi_end=360., dazi=1.,
179 phases={'P': 'first{stored:any_P}-10%',
180 'S': 'last{stored:any_S}+50'},
181 interpolation='multilinear',
182 target_depth=0.0,
183 quantity='displacement', envelope=False,
184 component='R', fmin=0.01, fmax=0.1,
185 hillshade=True, cmap=None,
186 plot_mt='full', show_phases=True, show_description=True,
187 reverse_time=False, show_nucleations=True, axes=None, nthreads=0):
188 '''
189 Plot the directivity and radiation characteristics of source models.
191 Synthetic seismic traces (R, T or Z) are forward-modelled at a defined
192 radius, covering the full or partial azimuthal range and projected on a
193 polar plot. Difference in the amplitude are enhanced by hillshading
194 the data.
196 :param engine: Forward modelling engine
197 :type engine: :py:class:`~pyrocko.gf.seismosizer.Engine`
198 :param source: Parametrized source model
199 :type source: :py:class:`~pyrocko.gf.seismosizer.Source`
200 :param store_id: Store ID used for forward modelling
201 :type store_id: str
202 :param distance: Distance in [m]
203 :type distance: float
204 :param azi_begin: Begin azimuth in [deg]
205 :type azi_begin: float
206 :param azi_end: End azimuth in [deg]
207 :type azi_end: float
208 :param dazi: Delta azimuth, bin size [deg]
209 :type dazi: float
210 :param phases: Phases to define start and end of time window
211 :type phases: :py:class:`dict` with :py:class:`str` keys and
212 :py:class:`~pyrocko.gf.meta.Timing` values
213 :param quantity: Seismogram quantity, default ``displacement``
214 :type quantity: str
215 :param envelope: Plot envelope 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: Matplotlib 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 description, 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 depth=target_depth,
264 interpolation=interpolation,
265 components='R', quantity=quantity)
266 ref_target = targets[0]
267 store = engine.get_store(store_id)
268 mt = source.pyrocko_moment_tensor(store=store, target=ref_target)
270 resp = engine.process(source, targets, nthreads=nthreads)
271 data, times = get_seismogram_array(
272 resp, fmin, fmax,
273 component=component, envelope=envelope)
275 nucl_depth = source.depth
276 nucl_distance = distance
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.warning('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 ax.grid(False)
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)