Browse Source

directivity plot: fixes and improvements

- fix ignored `quantity` param
- fix color scale range problem
- add `interpolation` param to set GF interpolation methd
- add `target_depth` param to set receiver depth
- changed waveform processing to prevent demeaning of traces
- changed waveform processing to not do bandpass but separate high and low
- docs: add simple analytical radiation pattern example
- fix typos
Sebastian Heimann 1 month ago
parent
commit
df1f097b8b
  1. 37
      doc/source/library/examples/plotting.rst
  2. BIN
      doc/source/static/directivity_envelope_rectangular.png
  3. BIN
      doc/source/static/directivity_rectangular.png
  4. BIN
      doc/source/static/radiation_pattern.png
  5. 19
      examples/plot_directivity.py
  6. 70
      examples/plot_radiation_pattern.py
  7. 8
      src/fomosto/ahfullgreen.py
  8. 30
      src/plot/directivity.py

37
doc/source/library/examples/plotting.rst

@ -50,9 +50,9 @@ Classes covered in these examples:
* :class:`pyrocko.gf.seismosizer.DCSource` (a representation of a double
couple source object),
* :class:`pyrocko.gf.seismosizer.RectangularExplosionSource` (a
representation of a rectangular explostion source),
representation of a rectangular explosion source),
* :class:`pyrocko.gf.seismosizer.CLVDSource` (a representation of a
compensated linear vector diploe source object)
compensated linear vector dipole source object)
* :class:`pyrocko.gf.seismosizer.DoubleDCSource` (a representation of a
double double-couple source object).
@ -96,7 +96,7 @@ Beachballs from source objects
This example shows how to add beachballs of various sizes to the corners of a
plot by obtaining the moment tensor from four different source object types:
:class:`pyrocko.gf.seismosizer.DCSource` (upper left),
:class:`pyrocko.gf.seismosizer.RectangularExplosionSource` (upper right),
:class:`pyrocko.gf.seismosizer.RectangularExplosionSource` (upper right),
:class:`pyrocko.gf.seismosizer.CLVDSource` (lower left) and
:class:`pyrocko.gf.seismosizer.DoubleDCSource` (lower right).
@ -186,7 +186,7 @@ components.
The function :py:func:`pyrocko.plot.hudson.project` may be used to get the
*(u,v)* coordinates for a given (full) moment tensor used for positioning the
symbol in the plot. The function :py:func:`pyrocko.plot.hudson.draw_axes` can
be used to conveniently draw the axes and annotions. Note, that we follow the
be used to conveniently draw the axes and annotations. Note, that we follow the
original convention introduced by Hudson, to place the negative CLVD on the
right hand side.
@ -205,10 +205,37 @@ Download :download:`hudson_diagram.py </../../examples/hudson_diagram.py>`
Source radiation plot
---------------------
The directivity and radiation characteristics of any point or finite
The directivity and radiation characteristics of any point or finite
:py:class:`~pyrocko.gf.seismosizer.Source` model can be illustrated with
:py:func:`~pyrocko.plot.directivity.plot_directivity`.
Radiation pattern effects
^^^^^^^^^^^^^^^^^^^^^^^^^
The following educational example illustrates radiation pattern effects from a
point source in a homogeneous full space. Analytical Green's functions for a
homogeneous full space are computed within the example code by use of the
ahfullgreen backend of Fomosto.
Download
:download:`plot_directivity.py </../../examples/plot_radiation_pattern.py>`
.. literalinclude :: /../../examples/plot_radiation_pattern.py
:language: python
.. figure :: /static/radiation_pattern.png
:align: center
:alt: Source radiation pattern of a double-couple point source in
a homogeneous full space.
Radial component radiation pattern for a double-couple point source in a
homogeneous full space observed in the plane of the source at a distance of
10 km. Note that the S waves seen in this example are pure near-field
effects. They get less pronounced when going to higher frequencies.
Directivity effects
^^^^^^^^^^^^^^^^^^^
Synthetic seismic traces (R, T or Z) are forward-modelled at a defined radius,
covering the full or partial azimuthal range and projected on a polar plot.
Difference in the amplitude are enhanced by hillshading the data.

BIN
doc/source/static/directivity_envelope_rectangular.png

Binary file not shown.

Before

Width:  |  Height:  |  Size: 237 KiB

After

Width:  |  Height:  |  Size: 184 KiB

BIN
doc/source/static/directivity_rectangular.png

Binary file not shown.

Before

Width:  |  Height:  |  Size: 243 KiB

After

Width:  |  Height:  |  Size: 191 KiB

BIN
doc/source/static/radiation_pattern.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 147 KiB

19
examples/plot_directivity.py

@ -24,22 +24,31 @@ rect_source = RectangularSource(
dip=76.6,
rake=-.4,
anchor='top',
nucleation_x=-.57,
nucleation_y=-.59,
velocity=2070.,
length=27*km,
width=9.4*km,
slip=1.4)
# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(7,5))
# axes = fig.add_subplot(111, polar=True)
resp = plot_directivity(
engine, rect_source, store_id,
distance=300*km, dazi=5., component='R',
plot_mt='full', show_phases=True,
# axes=axes,
distance=300*km,
dazi=5.,
component='R',
plot_mt='full',
show_phases=True,
phases={
'First': 'first{stored:begin}-10%',
'Last': 'last{stored:end}+20'
},
quantity='displacement', envelope=True)
quantity='displacement',
envelope=False)
# fig.savefig('directivity_rectangular.png')

70
examples/plot_radiation_pattern.py

@ -0,0 +1,70 @@
import os
import shutil
from pyrocko.plot.directivity import plot_directivity
from pyrocko.gf import LocalEngine, DCSource, Store
from pyrocko.fomosto import ahfullgreen
km = 1e3
def make_homogeneous_gf_store(
path, store_id, source_depth, receiver_depth, distance):
if os.path.exists(path):
shutil.rmtree(path)
ahfullgreen.init(path, None, config_params=dict(
id=store_id,
sample_rate=20.,
receiver_depth=receiver_depth,
source_depth_min=source_depth,
source_depth_max=source_depth,
distance_min=distance,
distance_max=distance))
store = Store(path)
store.make_ttt()
ahfullgreen.build(path)
store_id = 'gf_homogeneous_radpat'
store_path = os.path.join('.', store_id)
distance = 10*km
receiver_depth = 0.0
source = DCSource(
depth=0.,
strike=0.,
dip=90.,
rake=0.)
make_homogeneous_gf_store(
store_path, store_id, source.depth, receiver_depth, distance)
engine = LocalEngine(store_dirs=[store_path])
# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(7,5))
# axes = fig.add_subplot(111, polar=True)
resp = plot_directivity(
engine, source, store_id,
# axes=axes,
distance=distance,
dazi=5.,
component='R',
target_depth=receiver_depth,
plot_mt='full',
show_phases=True,
fmin=None,
fmax=1.0,
phases={
'P': '{stored:anyP}-50%',
'S': '{stored:anyS}+50%'
},
interpolation='nearest_neighbor',
quantity='velocity',
envelope=False,
hillshade=False)
# fig.savefig('radiation_pattern.png')

8
src/fomosto/ahfullgreen.py

@ -211,14 +211,14 @@ class AhfullGFBuilder(gf.builder.Builder):
(index+1, self.nblocks))
def init(store_dir, variant):
def init(store_dir, variant, config_params=None):
assert variant is None
modelling_code_id = 'ahfullgreen'
store_id = os.path.basename(os.path.realpath(store_dir))
config = gf.meta.ConfigTypeA(
d = dict(
id=store_id,
ncomponents=10,
sample_rate=20.,
@ -239,6 +239,10 @@ def init(store_dir, variant):
id='anyS',
definition='S,s,\\S,\\s')])
if config_params is not None:
d.update(config_params)
config = gf.meta.ConfigTypeA(**d)
config.validate()
return gf.store.Store.create_editables(
store_dir, config=config)

30
src/plot/directivity.py

@ -33,6 +33,7 @@ QUANTITY_LABEL = {
def get_azimuthal_targets(
store_id, source, radius,
azi_begin=0., azi_end=360., dazi=1.,
depth=0.0,
interpolation='multilinear',
components='RTZ', quantity='displacement'):
@ -54,7 +55,7 @@ def get_azimuthal_targets(
assert comp in dips.keys()
target_kwargs = dict(
quantity='displacement',
quantity=quantity,
interpolation=interpolation,
store_id=store_id)
@ -64,6 +65,7 @@ def get_azimuthal_targets(
lon=source.lon,
north_shift=coords[0, iazi] + source.north_shift,
east_shift=coords[1, iazi] + source.east_shift,
depth=depth,
azimuth={
'R': azi,
'T': azi+90.,
@ -98,12 +100,11 @@ def get_seismogram_array(
assert hasattr(target, 'azimuth')
assert target.dazi
if fmin and fmax:
tr.bandpass(2, fmin, fmax)
elif fmin:
tr.highpass(4, fmin)
elif fmax:
tr.lowpass(4, fmax)
if fmin is not None:
tr.highpass(4, fmin, demean=False)
if fmax is not None:
tr.lowpass(4, fmax, demean=False)
tmin = min(tmin, tr.tmin) if tmin else tr.tmin
tmax = max(tmax, tr.tmax) if tmax else tr.tmax
@ -115,7 +116,6 @@ def get_seismogram_array(
tr.abshilbert()
data = num.array([tr.get_ydata() for tr in traces])
data -= data.mean()
nsamples = data.shape[1]
return data, num.linspace(tmin, tmax, nsamples)
@ -179,6 +179,8 @@ def plot_directivity(
distance=300*km, azi_begin=0., azi_end=360., dazi=1.,
phases={'P': 'first{stored:any_P}-10%',
'S': 'last{stored:any_S}+50'},
interpolation='multilinear',
target_depth=0.0,
quantity='displacement', envelope=False,
component='R', fmin=0.01, fmax=0.1,
hillshade=True, cmap=None,
@ -210,9 +212,11 @@ def plot_directivity(
:type phase_begin: :py:class:`~pyrocko.gf.meta.Timing`
:param phase_end: End time of the window
:type phase_end: :py:class:`~pyrocko.gf.meta.Timing`
:param target_depth: Receiver depth [m]
:type target_depth: float
:param quantity: Seismogram quantity, default ``displacement``
:type quantity: str
:param envelope: Plot envelop instead of seismic trace
:param envelope: Plot envelope instead of seismic trace
:type envelope: bool
:param component: Forward modelled component, default ``R``. Choose from
`RTZ`
@ -223,7 +227,7 @@ def plot_directivity(
:type fmax: float
:param hillshade: Enable hillshading, default ``True``
:type hillshade: bool
:param cmap: Matplotlit colormap to use, default ``seismic``.
:param cmap: Matplotlib colormap to use, default ``seismic``.
When ``envelope`` is ``True`` the default colormap will be ``Reds``.
:type cmap: str
:param plot_mt: Plot a centered moment tensor, default ``full``.
@ -231,7 +235,7 @@ def plot_directivity(
:type plot_mt: str, bool
:param show_phases: Show annotations, default ``True``
:type show_phases: bool
:param show_description: Show desciption, default ``True``
:param show_description: Show description, default ``True``
:type show_description: bool
:param reverse_time: Reverse time axis. First phases arrive at the center,
default ``False``
@ -260,6 +264,8 @@ def plot_directivity(
targets, azimuths = get_azimuthal_targets(
store_id, source, distance, azi_begin, azi_end, dazi,
depth=target_depth,
interpolation=interpolation,
components='R', quantity=quantity)
ref_target = targets[0]
store = engine.get_store(store_id)
@ -366,7 +372,7 @@ def plot_directivity(
mesh = ax.pcolormesh(
azimuths * d2r, times, data,
cmap=cmw.cmap, shading='gouraud', zorder=0)
cmap=cmw.cmap, norm=cmw.norm, shading='gouraud', zorder=0)
if hillshade:
mesh.update_scalarmappable()

Loading…
Cancel
Save