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

CHANGELOG.md

@@ -32,6 +32,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - Fixed `Trace.envelope`.
 
 ### Changed
+- Improvements to directivity plot.
 - Renamed Snuffling "Block RMS" to "RMS".
 - Fomosto QSEIS backend can now handle non-zero (negative) top layer depth,
   e.g. for models including the atmosphere.

doc/source/library/examples/plotting.rst

@@ -65,9 +65,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).
 
@@ -111,7 +111,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).
 
@@ -201,7 +201,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.
 
@@ -220,10 +220,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.

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')

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')

src/fomosto/ahfullgreen.py

@@ -220,7 +220,7 @@ 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'
@@ -229,7 +229,7 @@ def init(store_dir, variant):
 
     ahfullgreen = AhfullgreenConfig()
 
-    config = gf.meta.ConfigTypeA(
+    d = dict(
         id=store_id,
         ncomponents=10,
         sample_rate=20.,
@@ -250,6 +250,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, extra={'ahfullgreen': ahfullgreen})

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.envelope()
 
     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,
@@ -211,7 +213,7 @@ def plot_directivity(
         :py:class:`~pyrocko.gf.meta.Timing` values
     :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`
@@ -222,7 +224,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``.
@@ -230,7 +232,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``
@@ -259,6 +261,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)
@@ -272,8 +276,6 @@ def plot_directivity(
     nucl_depth = source.depth
     nucl_distance = distance
 
-    anch_x, anch_y = map_anchor[source.anchor]
-
     if hasattr(source, 'nucleation_x') and hasattr(source, 'nucleation_y'):
         try:
             iter(source.nucleation_x)
@@ -367,7 +369,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()