add okada modelling and pseudo-dynamic rupture model (chapter 11)

- docs
- examples
- polishing
- bugfixes

examples/automap_example2.py

@@ -55,6 +55,7 @@ def extend_1d_coordinate_array(array):
 
     return out
 
+
 def extend_2d_data_array(array):
     '''
     Extend 2D data array for gridded data, that grid corners are plotted
@@ -72,6 +73,7 @@ def extend_2d_data_array(array):
 
     return out
 
+
 def tile_to_length_width(m, ref_lat, ref_lon):
     '''
     Transform grid tile (lat, lon) to easting, northing for data interpolation
@@ -88,6 +90,7 @@ def tile_to_length_width(m, ref_lat, ref_lon):
     return num.hstack((
         grid_easting.reshape(-1, 1), grid_northing.reshape(-1, 1)))
 
+
 def data_to_grid(m, x, y, data):
     '''
     Create data grid from data and coordinate arrays

examples/gf_forward_pseudo_rupture.py

@@ -44,9 +44,9 @@ source_params = dict(
 dyn_rupture = gf.PseudoDynamicRupture(
     nx=5, ny=5,
     tractions=gf.tractions.HomogeneousTractions(
-        t_strike=1.e4,
-        t_dip=0.e4,
-        t_normal=0.),
+        strike=1.e4,
+        dip=0.e4,
+        normal=0.),
     **source_params)
 
 dyn_rupture.discretize_patches(store)

examples/okada_forward_example.py

@@ -6,8 +6,8 @@ from pyrocko.plot import dislocation as displt
 d2r = num.pi / 180.
 km = 1000.
 
-# Set source parameter
-src_north, src_east, src_depth = 0. * km, 0. * km, 100. * km
+# Set source parameters
+src_north, src_east, src_depth = 20. * km, -45. * km, 10. * km
 
 length_total = 50. * km
 width_total = 15. * km
@@ -24,7 +24,7 @@ source = OkadaSource(
     lat=0., lon=0., north_shift=src_north, east_shift=src_east,
     depth=src_depth,
     al1=al1, al2=al2, aw1=aw1, aw2=aw2,
-    strike=45., dip=90., rake=90.,
+    strike=66., dip=45., rake=90.,
     slip=1., opening=0., poisson=0.25, shearmod=32.0e9)
 
 source_discretized, _ = source.discretize(nlength, nwidth)
@@ -44,7 +44,8 @@ source_disl = num.array([
     patch.source_disloc() for patch in source_discretized])
 result = okada_ext.okada(
     source_patch, source_disl, receiver_coords,
-    source.lamb, source.shearmod, 0)
+    source.lamb, source.shearmod, nthreads=0, rotate_sdn=False,
+    stack_sources=True)
 
 # Plot
 displt.plot(result, receiver_coords, cmap='coolwarm', zero_center=True)

examples/okada_inversion_example.py

@@ -34,7 +34,7 @@ source_discretized, _ = source.discretize(nlength, nwidth)
 receiver_coords = num.array([
     src.source_patch()[:3] for src in source_discretized])
 
-# Create Stress drop (traction) array
+# Create stress drop (traction) array with spatial varying traction vectors
 dstress = -1.5e6
 stress_comp = 1
 
@@ -50,10 +50,11 @@ for il in range(nlength):
                 (iw > 2 and iw < nwidth - 4):
             stress_field[idx + stress_comp] = dstress / 4.
 
+# Invert for dislocation on source plane based on given stress field
 disloc_est = DislocationInverter.get_disloc_lsq(
-    stress_field, source_list=source_discretized)
-
+    stress_field, source_list=source_discretized, nthreads=0)
 
+# Plot
 displt.plot(
     disloc_est.reshape(npoints, 3),
     receiver_coords,

src/ahfullgreen.py

@@ -7,9 +7,13 @@ from __future__ import absolute_import, division
 import math
 import numpy as num
 from . import trace
+from .guts import Float, Object
 from . import ahfullgreen_ext as ext
 
 
+guts_prefix = 'pf'
+
+
 class AhfullgreenError(Exception):
     pass
 
@@ -21,7 +25,7 @@ def make_seismogram(
         npad_levelling=40, out_alignment=0.):
 
     if stf is None:
-        stf = Impulse()
+        stf = AhfullgreenSTFImpulse()
 
     x = num.asarray(x, float)
     f = num.asarray(f, float)
@@ -157,10 +161,11 @@ def add_seismogram(
         out[ind(tmax, out_offset)+1:] += temp[ind(tmax, tstart)]
 
 
-class Impulse(object):
+class AhfullgreenSTF(Object):
+    pass
 
-    def __init__(self):
-        pass
+
+class AhfullgreenSTFImpulse(AhfullgreenSTF):
 
     def t_cutoff(self):
         return None
@@ -169,14 +174,14 @@ class Impulse(object):
         return num.ones(f.size, dtype=complex)
 
 
-class Gauss(object):
-    def __init__(self, tau):
-        self._tau = tau
+class AhfullgreenSTFGauss(AhfullgreenSTF):
+
+    tau = Float.T(default=1.0)
 
     def t_cutoff(self):
-        return self._tau * 2.
+        return self.tau * 2.
 
     def __call__(self, f):
         omega = f * 2. * math.pi
 
-        return num.exp(-(omega**2 * self._tau**2 / 8.))
+        return num.exp(-(omega**2 * self.tau**2 / 8.))

src/fomosto/ahfullgreen.py

@@ -10,8 +10,10 @@ import os
 import math
 import signal
 
+from pyrocko.guts import Object
 from pyrocko import trace, cake, gf
-from pyrocko.ahfullgreen import add_seismogram, Impulse, Gauss
+from pyrocko.ahfullgreen import add_seismogram, AhfullgreenSTFImpulse, \
+    AhfullgreenSTF
 from pyrocko.moment_tensor import MomentTensor, symmat6
 
 km = 1000.
@@ -41,8 +43,12 @@ class AhfullgreenError(gf.store.StoreError):
     pass
 
 
+class AhfullgreenConfig(Object):
+    stf = AhfullgreenSTF.T(default=AhfullgreenSTFImpulse.D())
+
+
 def make_traces(material, source_mech, deltat, norths, easts,
-                source_depth, receiver_depth):
+                source_depth, receiver_depth, stf):
 
     if isinstance(source_mech, MomentTensor):
         m6 = source_mech.m6()
@@ -72,7 +78,7 @@ def make_traces(material, source_mech, deltat, norths, easts,
             material.vp, material.vs, material.rho, material.qp, material.qs,
             x, f, m6, 'displacement',
             deltat, tmin, outx, outy, outz,
-            stf=Gauss(tau=1.))
+            stf=stf)
 
         for i_comp, o in enumerate((outx, outy, outz)):
             comp = components[i_comp]
@@ -101,6 +107,8 @@ class AhfullGFBuilder(gf.builder.Builder):
         gf.builder.Builder.__init__(
             self, self.store.config, step, block_size=block_size, force=force)
 
+        self.ahfullgreen_config = self.store.get_extra('ahfullgreen')
+
     def cleanup(self):
         self.store.close()
 
@@ -158,7 +166,8 @@ class AhfullGFBuilder(gf.builder.Builder):
             rawtraces = make_traces(
                 self.store.config.earthmodel_1d.require_homogeneous(),
                 source_mech, 1.0/self.store.config.sample_rate,
-                distances, num.zeros_like(distances), sz, rz)
+                distances, num.zeros_like(distances), sz, rz,
+                self.ahfullgreen_config.stf)
 
             interrupted = []
 
@@ -218,6 +227,8 @@ def init(store_dir, variant):
 
     store_id = os.path.basename(os.path.realpath(store_dir))
 
+    ahfullgreen = AhfullgreenConfig()
+
     config = gf.meta.ConfigTypeA(
         id=store_id,
         ncomponents=10,
@@ -241,7 +252,7 @@ def init(store_dir, variant):
 
     config.validate()
     return gf.store.Store.create_editables(
-        store_dir, config=config)
+        store_dir, config=config, extra={'ahfullgreen': ahfullgreen})
 
 
 def build(store_dir, force=False, nworkers=None, continue_=False, step=None,

src/modelling/ext/okada_ext.c

@@ -853,10 +853,37 @@ static okada_error_t dc3d_flexi(
         int rot_sdn) {
 
     /*
-     * Wrapper for dc3d function
-     * Is applied on single source-receiver pair
-     * source is arbitrary oriented
-    */
+     *   Wrapper for dc3d function
+     *   It is applied on single source-receiver pair for an arbitrary oriented source
+     *   Displacement and strain at depth due to buried finite fault in a
+     *    semiinfinite medium.
+     *
+     *    Input:
+     *
+     *      alpha:       medium constant  (lambda + myu)/(lambda + 2*myu)
+     *      nr, er, zr:  coordinate of observing point north, east, down [m]
+     *      ns, es, zs:  coordinate of source plane reference point north, east, down [m]
+     *      strike:      strike-angle (degree)
+     *      dip:         dip-angle (degree)
+     *      al1, al2:    fault length range
+     *      aw1, aw2:    fault width range
+     *      disl1-disl3: strike-, dip-, tensile-dislocations
+     *      rot_sdn:     1, if displacements and derivatives are returned in strike-dip-
+     *                   tensile (normal), 0 for north-east-down
+     *    Output:
+     *
+     *      u[12]: displacement (units of disl) and derivatives
+     *             ((unit of disl) / (unit of nr, er, zr, ns, es, zs, al, aw)) as
+     *
+     *           [us, ud, un, uss, uds, uns, usd, udd, und, usn, udn, unn]
+     *              with rot_sdn = 1 and s-strike, d-down dip and n-normal
+     *           [un, ue, ud, unn, uen, udn, une, uee, ude, und, ued, udd]
+     *              with rot_sdn = 0 and n-north, e-east and d-down
+     *
+     *    Return value:
+     *
+     *      0 for success or error code
+     */
 
     double rotmat[3][3];
     double r[3], rrot[3];
@@ -874,11 +901,15 @@ static okada_error_t dc3d_flexi(
 
     rot_vec31(r, rotmat, rrot);
 
+    if (dr == 0.){
+        rrot[2] = 0.;
+    }
+
     if (dip == 90.) {
         dip -= 1E-2;
     }
     iret = dc3d(alpha, rrot[0], rrot[1], rrot[2], ds, dip, al1, al2, aw1, aw2, disl1, disl2, disl3, uokada);
-    
+
     /*
      * Back rotation of displacement and strain vector/tensor
      * with the transposed rotmat equals rotmat

src/modelling/okada.py

@@ -314,7 +314,7 @@ class DislocationInverter(object):
             or additionally include opening, default False.
         :type pure_shear: optional, bool
         :param rotate_sdn: Rotation towards strike, dip, normal, default True.
-        :type pure_shear: optional, bool
+        :type rotate_sdn: optional, bool
         :param nthreads: Number of threads, default 1
         :type nthreads: optional, int
 
@@ -410,7 +410,11 @@ class DislocationInverter(object):
             :py:class:`pyrocko.modelling.OkadaSource`
         :param pure_shear: Flag, if also opening mode shall be taken into
             account (False) or the fault is described as pure shear (True).
-        :type pure_shear: optional, Bool
+        :type pure_shear: optional, bool
+        :param rotate_sdn: Rotation towards strike, dip, normal, default True.
+        :type rotate_sdn: optional, bool
+        :param nthreads: Number of threads, default 1
+        :type nthreads: optional, int
 
         :return: coefficient matrix for all sources
         :rtype: :py:class:`numpy.ndarray`,
@@ -505,7 +509,11 @@ class DislocationInverter(object):
             :py:class:`pyrocko.modelling.OkadaSource`
         :param pure_shear: Flag, if also opening mode shall be taken into
             account (False) or the fault is described as pure shear (True).
-        :type pure_shear: optional, Bool
+        :type pure_shear: optional, bool
+        :param rotate_sdn: Rotation towards strike, dip, normal, default True.
+        :type rotate_sdn: optional, bool
+        :param nthreads: Number of threads, default 1
+        :type nthreads: optional, int
 
         :return: coefficient matrix for all sources
         :rtype: :py:class:`numpy.ndarray`,

src/plot/dislocation.py

@@ -41,8 +41,9 @@ def draw(
             num.min(dislocation), num.max(dislocation)]))
         vmin = -vmax
     else:
-        vmax = num.max(dislocation)
         vmin = num.min(dislocation)
+        vmax = num.max(dislocation)
+
     scat = axes.scatter(
         coordinates[:, 1],
         coordinates[:, 0],

src/plot/dynamic_rupture.py

@@ -16,27 +16,29 @@ from scipy.interpolate import RegularGridInterpolator as scrgi
 
 from matplotlib import cm, pyplot as plt
 
-from pyrocko import gmtpy, moment_tensor as pmt, orthodrome as pod, util
+from pyrocko import gmtpy, orthodrome as pod
 from pyrocko.plot import (mpl_init, mpl_margins, mpl_papersize, mpl_color,
                           AutoScaler)
-from pyrocko.plot.automap import Map
+from pyrocko.plot.automap import Map, NoTopo
 from pyrocko.gf import PseudoDynamicRupture
 from pyrocko.gf.seismosizer import map_anchor
+from pyrocko.dataset.topo.tile import Tile
 
 logger = logging.getLogger('pyrocko.plot.dynamic_rupture')
 
+gmtpy.check_have_gmt()
+gmt = gmtpy.GMT()
 
 km = 1e3
 
 d2m = 111180.
 m2d = 1. / d2m
 
+cm2inch = gmtpy.cm/gmtpy.inch
+
 d2r = num.pi / 180.
 r2d = 1. / d2r
 
-gmtpy.check_have_gmt()
-gmt = gmtpy.GMT()
-
 
 def _save_grid(lats, lons, data, filename):
     '''
@@ -55,7 +57,7 @@ def _save_grid(lats, lons, data, filename):
     gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat')
 
 
-def _mplcmap_to_gmtcpt_code(mplcmap):
+def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
     '''
     Get gmt readable R/G/A code from a given matplotlib colormap
 
@@ -68,7 +70,7 @@ def _mplcmap_to_gmtcpt_code(mplcmap):
 
     cmap = cm.get_cmap(mplcmap)
 
-    rgbas = [cmap(i) for i in num.linspace(0, 255, 256).astype(num.int64)]
+    rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)]
 
     return ','.join(['%g/%g/%g' % (
         c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas])
@@ -104,6 +106,9 @@ def make_colormap(
     incr = scale[2]
     margin = 0.
 
+    if vmin == vmax:
+        space = True
+
     if space:
         margin = incr
 
@@ -125,7 +130,7 @@ def make_colormap(
 
     return gmt.makecpt(
         C=C,
-        D='i',
+        D='o',
         T='%g/%g/%g' % (
             vmin - margin, vmax + margin, incr),
         Z=True,
@@ -175,11 +180,13 @@ def xy_to_latlon(source, x, y):
     s = source
     ax, ay = map_anchor[s.anchor]
 
-    l, w = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
+    length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
     strike, dip = s.strike * d2r, s.dip * d2r
 
-    northing = num.cos(strike) * l - num.cos(dip) * num.sin(strike) * w
-    easting = num.sin(strike) * l + num.cos(dip) * num.cos(strike) * w
+    northing = \
+        num.cos(strike) * length - num.cos(dip) * num.sin(strike) * width
+    easting = \
+        num.sin(strike) * length + num.cos(dip) * num.cos(strike) * width
 
     return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
 
@@ -200,12 +207,12 @@ def xy_to_lw(source, x, y):
     :rtype: tuple of float
     '''
 
-    l, w = source.length, source.width
+    length, width = source.length, source.width
 
     ax, ay = map_anchor[source.anchor]
 
-    lengths = (x - ax) / 2. * l
-    widths = (y - ay) / 2. * w
+    lengths = (x - ax) / 2. * length
+    widths = (y - ay) / 2. * width
 
     return lengths, widths
 
@@ -240,18 +247,18 @@ cbar_helper = {
         'factor': 1.},
     'strike': {
         'unit': '°',
-        'factor': 1e-6},
+        'factor': 1.},
     'dip': {
         'unit': '°',
-        'factor': 1e-6},
+        'factor': 1.},
     'vr': {
         'unit': 'km/s',
         'factor': 1e-3},
     'length': {
-        'unit': 'm',
+        'unit': 'km',
         'factor': 1e-3},
     'width': {
-        'unit': 'm',
+        'unit': 'km',
         'factor': 1e-3}
 }
 
@@ -314,21 +321,28 @@ class RuptureMap(Map):
             fontcolor='darkslategrey',
             width=20.,
             height=14.,
+            margins=None,
+            topo_cpt_wet='light_sea_uniform',
+            topo_cpt_dry='light_land_uniform',
             *args, **kwargs):
 
         size = (width, height)
         fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
-        margins = [
-            fontsize * 0.15, num.min(size) / 200.,
-            num.min(size) / 200., fontsize * 0.05]
+
+        if margins is None:
+            margins = [
+                fontsize * 0.15, num.min(size) / 200.,
+                num.min(size) / 200., fontsize * 0.05]
 
         Map.__init__(self, margins=margins, width=width, height=height,
                      gmt_config=gmt_config,
+                     topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
                      *args, **kwargs)
 
         self.source = source
         self._fontcolor = fontcolor
         self._fontsize = fontsize
+        self.axes_layout = 'WeSn'
 
     @property
     def size(self):
@@ -376,7 +390,22 @@ class RuptureMap(Map):
         if self._dems is None:
             self._setup()
 
-        t, _ = self._get_topo_tile('land')
+        try:
+            t, _ = self._get_topo_tile('land')
+        except NoTopo:
+            wesn = self.wesn
+
+            nx = int(num.ceil(
+                self.width * cm2inch * self.topo_resolution_max))
+            ny = int(num.ceil(
+                self.height * cm2inch * self.topo_resolution_max))
+
+            data = num.zeros((nx, ny))
+
+            t = Tile(wesn[0], wesn[2],
+                     (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
+                     data)
+
         return t
 
     def _patches_to_lw(self):
@@ -399,7 +428,7 @@ class RuptureMap(Map):
         patch_lengths = num.concatenate((
             num.array([0.]),
             num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
-            num.array([patch_l* source.nx])))
+            num.array([patch_l * source.nx])))
 
         patch_widths = num.concatenate((
             num.array([0.]),
@@ -719,6 +748,33 @@ class RuptureMap(Map):
             *self.jxyr,
             **kwargs)
 
+    def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
+        ''' Draw vectors based on two grid files
+
+        Two grid files for vector lengths in x and y need to be given. The
+        function calls gmt.grdvector. All arguments defined for this function
+        in gmt can be passed as keyword arguments. Different standard settings
+        are applied if not defined differently.
+
+        :param x_gridfile: File of the grid defining vector lengths in x
+        :type x_gridfile: string
+        :param y_gridfile: File of the grid defining vector lengths in y
+        :type y_gridfile: string
+        :param vcolor: Vector face color as defined in "G" option
+        :type vcolor: string
+        '''
+
+        kwargs['S'] = kwargs.get('S', 'il1.')
+        kwargs['I'] = kwargs.get('I', 'x20')
+        kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
+        kwargs['Q'] = kwargs.get('Q', '4c+e+n5c+h1')
+
+        self.gmt.grdvector(
+            x_gridfile, y_gridfile,
+            G='%s' % 'lightgrey' if not vcolor else vcolor,
+            *self.jxyr,
+            **kwargs)
+
     def draw_dynamic_data(self, data, **kwargs):
         '''
         Draw an image of any data gridded on the patches e.g dislocation
@@ -761,17 +817,13 @@ class RuptureMap(Map):
         source = self.source
 
         if a == 'traction':
-            data = num.linalg.norm(source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny), axis=1)
+            data = num.linalg.norm(source.get_tractions(), axis=1)
         elif a == 'tx':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 0]
+            data = source.get_tractions()[:, 0]
         elif a == 'ty':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 1]
+            data = source.get_tractions()[:, 1]
         elif a == 'tz':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 2]
+            data = source.get_tractions()[:, 2]
         else:
             data = source.get_patch_attribute(attribute)
 
@@ -841,15 +893,8 @@ class RuptureMap(Map):
             point='p',
             triangle='t')
 
-        try:
-            iterator = iter(lats)
-        except TypeError:
-            lats = num.array([lats])
-
-        try:
-            iterator = iter(lons)
-        except TypeError:
-            lons = num.array([lons])
+        lats = num.atleast_1d(lats)
+        lons = num.atleast_1d(lons)
 
         if lats.shape[0] != lons.shape[0]:
             raise IndexError('lats and lons do not have the same shape!')
@@ -937,6 +982,36 @@ class RuptureMap(Map):
 
         clear_temp(gridfiles=[tmp_grd_file], cpts=[])
 
+    def draw_dislocation_vector(self, time=None, **kwargs):
+        '''Draw vector arrows onto map indicating direction of dislocation
+
+        For a given time (if time is None, tmax is used) and given component
+        the dislocation is plotted as vectors onto the map.
+        :param time: time after origin, for which dislocation is computed. If
+            None, tmax is taken.
+        :type time: optional, float
+        '''
+
+        disl = self.source.get_okada_slip(time=time)
+
+        p_strike = self.source.get_patch_attribute('strike') * d2r
+        p_dip = self.source.get_patch_attribute('dip') * d2r
+
+        disl_dh = num.cos(p_dip) * disl[:, 1]
+        disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
+        disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
+
+        tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
+
+        self.patch_data_to_grid(disl_n, tmp_grd_files[0])
+        self.patch_data_to_grid(disl_e, tmp_grd_files[1])
+
+        self.draw_vector(
+            tmp_grd_files[1], tmp_grd_files[0],
+            **kwargs)
+
+        clear_temp(gridfiles=tmp_grd_files, cpts=[])
+
 
 class RuptureView(object):
     _patches_to_lw = RuptureMap._patches_to_lw
@@ -1035,12 +1110,12 @@ class RuptureView(object):
                 *args,
                 **kwargs)
 
-    def _draw_image(self, l, w, data, *args, **kwargs):
+    def _draw_image(self, length, width, data, *args, **kwargs):
         if self._axes is not None:
             if 'extent' not in kwargs:
                 kwargs['extent'] = [
-                    num.min(l), num.max(l),
-                    num.max(w), num.min(w)]
+                    num.min(length), num.max(length),
+                    num.max(width), num.min(width)]
 
             im = self._axes.imshow(
                 data,
@@ -1091,23 +1166,23 @@ class RuptureView(object):
                 *args,
                 **kwargs)
 
-    def draw_point(self, l, w, *args, **kwargs):
+    def draw_point(self, length, width, *args, **kwargs):
         ''' Draw a point onto the figure.
 
         Args and kwargs can be defined according to matplotlib.pyplot.scatter
-        :param l: Point(s) coordinate on the rupture plane along strike
+        :param length: Point(s) coordinate on the rupture plane along strike
             relative to the anchor point [m]
-        :type l: float, :py:class:`numpy.ndarray`
-        :param w: Point(s) coordinate on the rupture plane along downdip
+        :type length: float, :py:class:`numpy.ndarray`
+        :param width: Point(s) coordinate on the rupture plane along downdip
             relative to the anchor point [m]
-        :type w: float, :py:class:`numpy.ndarray`
+        :type width: float, :py:class:`numpy.ndarray`
         '''
 
         if self._axes is not None:
             kwargs['c'] = kwargs.get('c', mpl_color('scarletred2'))
             kwargs['s'] = kwargs.get('s', 100.)
 
-            self._axes.scatter(l, w, *args, **kwargs)
+            self._axes.scatter(length, width, *args, **kwargs)
 
     def draw_dynamic_data(self, data, **kwargs):
         ''' Draw an image of any data gridded on the patches e.g dislocation
@@ -1118,11 +1193,11 @@ class RuptureView(object):
 
         anchor_x, anchor_y = map_anchor[self.source.anchor]
 
-        l, w = xy_to_lw(
+        length, width = xy_to_lw(
             self.source, num.array([-1., 1.]), num.array([-1., 1.]))
 
-        l /= km
-        w /= km
+        length /= km
+        width /= km
 
         data = data.reshape(self.source.ny, self.source.nx, order='F')
 
@@ -1138,7 +1213,7 @@ class RuptureView(object):
                        ['xlabel', 'ylabel', 'title']])
 
         self._setup(**setup_kwargs)
-        self._draw_image(l=l, w=w, data=data, **kwargs)
+        self._draw_image(length=length, width=width, data=data, **kwargs)
 
     def draw_patch_parameter(self, attribute, **kwargs):
         ''' Draw an image of a chosen patch attribute e.g traction
@@ -1155,17 +1230,13 @@ class RuptureView(object):
         source = self.source
 
         if a == 'traction':
-            data = num.linalg.norm(source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny), axis=1)
+            data = num.linalg.norm(source.get_tractions(), axis=1)
         elif a == 'tx':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 0]
+            data = source.get_tractions()[:, 0]
         elif a == 'ty':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 1]
+            data = source.get_tractions()[:, 1]
         elif a == 'tz':
-            data = source.tractions.get_tractions(
-                nx=source.nx, ny=source.ny)[:, 2]
+            data = source.get_tractions()[:, 2]
         else:
             data = source.get_patch_attribute(attribute)
 
@@ -1196,15 +1267,15 @@ class RuptureView(object):
         if not clevel:
             clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
 
-        l, w = xy_to_lw(self.source, points_xy[:, 0], points_xy[:, 1])
+        length, width = xy_to_lw(self.source, points_xy[:, 0], points_xy[:, 1])
 
-        l /= km
-        w /= km
+        length /= km
+        width /= km
 
         kwargs['colors'] = kwargs.get('colors', '#474747ff')
         self._setup(**kwargs)
 
-        self._draw_contour(num.unique(l), num.unique(w), data=times.T,
+        self._draw_contour(num.unique(length), num.unique(width), data=times.T,
                            clevel=clevel, **kwargs)
 
     def draw_nucleation_point(self, **kwargs):
@@ -1213,12 +1284,12 @@ class RuptureView(object):
 
         nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
 
-        l, w = xy_to_lw(self.source, nuc_x, nuc_y)
+        length, width = xy_to_lw(self.source, nuc_x, nuc_y)
 
-        l /= km
-        w /= km
+        length /= km
+        width /= km
 
-        self.draw_point(l, w, marker='o', **kwargs)
+        self.draw_point(length, width, marker='o', **kwargs)
 
     def draw_dislocation(self, time=None, component='', **kwargs):
         ''' Draw dislocation onto map at any time
@@ -1276,16 +1347,16 @@ class RuptureView(object):
 
         anchor_x, anchor_y = map_anchor[self.source.anchor]
 
-        l, w = self._patches_to_lw()
-        l, w = l[1:-1], w[1:-1]
+        length, width = self._patches_to_lw()
+        length, width = length[1:-1], width[1:-1]
 
-        l /= km
-        w /= km
+        length /= km
+        width /= km
 
         kwargs['colors'] = kwargs.get('colors', '#474747ff')
 
         self._setup(**kwargs)
-        self._draw_contour(l, w, data=data, clevel=clevel, **kwargs)
+        self._draw_contour(length, width, data=data, clevel=clevel, **kwargs)
 
     def draw_source_dynamics(
             self, variable, store=None, dt=None, *args, **kwargs):
@@ -1311,7 +1382,7 @@ class RuptureView(object):
 
         v = variable
 
-        data, times = self.source.get_source_moment_rate(store=store, dt=dt)
+        data, times = self.source.get_moment_rate(store=store, deltat=dt)
 
         if v in ('moment_rate', 'stf'):
             name, unit = 'dM/dt', 'Nm/s'
@@ -1357,15 +1428,15 @@ class RuptureView(object):
         m = re.match(r'dislocation_([xyz])', v)
 
         if v in ('moment_rate', 'cumulative_moment', 'moment'):
-            data, times = source.get_moment_rate(dt=dt)
+            data, times = source.get_moment_rate_patches(dt=dt)
         elif 'dislocation' in v or 'slip_rate' == v:
             ddisloc, times = source.get_delta_slip(dt=dt)
 
         if v == 'moment_rate':
-            data, times = source.get_moment_rate(store=store, dt=dt)
+            data, times = source.get_moment_rate_patches(store=store, dt=dt)
             name, unit = 'dM/dt', 'Nm/s'
         elif v == 'cumulative_moment' or v == 'moment':
-            data, times = source.get_moment_rate(store=store, dt=dt)
+            data, times = source.get_moment_rate_patches(store=store, dt=dt)
             data = num.cumsum(data, axis=1)
             name, unit = 'M', 'Nm'
         elif v == 'slip_rate':
@@ -1454,6 +1525,38 @@ def render_movie(fn_path, output_path, dt=0.5):
         output_path])
 
 
+def render_gif(fn, output_path, loops=-1):
+    ''' Generate a gif based on a given mp4 using ffmpeg
+
+    Render a gif based on a given .mp4 movie file in fn_path.
+
+    :param fn: Path and file name of the input .mp4 file.
+    :type fn: string
+    :param output_path: Path and filename of the output gif movie file
+    :type output_path: string
+    :param loops: Number of gif repeats (loops). -1 is not repetition,
+        0 infinite
+    :type loops: optional, integer
+    '''
+
+    try:
+        check_call(['ffmpeg', '-loglevel', 'panic'])
+    except CalledProcessError:
+        pass
+    except (TypeError):
+        logger.warn(
+            'Package ffmpeg needed for movie rendering. Please install it '
+            '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
+        return
+
+    check_call([
+        'ffmpeg', '-i',
+        fn,
+        '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
+        '-loop', '%d' % loops,
+        output_path])
+
+
 def rupture_movie(
         source,
         store,
@@ -1463,10 +1566,13 @@ def rupture_movie(
         plot_type='map',
         dt=None,
         store_images=False,
+        render_as_gif=False,
+        gif_loops=-1,
+        original_patch_indices=None,
         **kwargs):
     ''' Generate a movie based on a given source for dynamic parameter
 
-    Create a .mp4 movie of one of the following dynamic parameters
+    Create a .mp4 movie or gif of one of the following dynamic parameters
     (dislocation, dislocation_x (along strike), dislocation_y (along updip),
     dislocation_z (normal), slip_rate, moment_rate). If whished, the single
     snap shots can be stored as images as well. kwargs have to be given
@@ -1486,7 +1592,7 @@ def rupture_movie(
     :type fn_path: optional, string
     :param prefix: File prefix used for the movie (and optional image) files
     :type prefix: optional, string
-    :param plot_type: Choice of plot type (map plot using
+    :param plot_type: Choice of plot type: 'map', 'view' (map plot using
         :py:class:`pyrocko.plot.dynamic_rupture.RuptureMap or plane view using
         :py:class:`pyrocko.plot.dynamic_rupture.RuptureView)
     :type plot_type: optional, string
@@ -1496,10 +1602,28 @@ def rupture_movie(
     :param store_images: Choice to store the single .png parameter snapshots in
         fn_path or not.
     :type store_images: optional, bool
+    :param render_as_gif: If True, the movie is converted into a gif. If False,
+        the movie is returned as mp4
+    :type render_as_gif: optional, bool
+    :param gif_loops: If render_as_gif is True, a gif with gif_loops number of
+        loops (repetitions) is returned. -1 is no repetition, 0 infinite.
+    :type gif_loops: optional, integer
+    :param original_patch_indices: If the source has been recalculated and
+        manually manipulated (patches, coef_mat etc.), then give here a boolean
+        array indicating the position of the remaining patches on the original
+        patch grid. False is interpreted as no data, True is replaced by the
+        patch data.
+    :type original_patch_indices: optional,
+        :py:class:`numpy.ndarray` of bool ``(source.nx, source.ny)``
     '''
 
     v = variable
-    source.discretize_basesource(store)
+
+    if not source.patches:
+        source.discretize_patches(store, interpolation='multilinear')
+
+    if source.coef_mat is None:
+        source.calc_coef_mat()
 
     if not prefix:
         prefix = v
@@ -1508,7 +1632,7 @@ def rupture_movie(
         dt = store.config.deltat
 
     if v == 'moment_rate':
-        data, times = source.get_moment_rate(dt=dt)
+        data, times = source.get_moment_rate_patches(dt=dt)
         name, unit = 'dM/dt', 'Nm/s'
     elif 'dislocation' in v or 'slip_rate' == v:
         ddisloc, times = source.get_delta_slip(dt=dt)
@@ -1519,10 +1643,10 @@ def rupture_movie(
 
     m = re.match(r'dislocation_([xyz])', v)
     if m:
-        data = num.cumsum(ddisloc, axis=2)[:, c2disl[m.group(1)], :]
+        data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
         name, unit = 'du%s' % m.group(1), 'm'
     elif v == 'dislocation':
-        data = num.linalg.norm(num.cumsum(ddisloc, axis=2), axis=1)
+        data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
         name, unit = 'du', 'm'
     elif v == 'slip_rate':
         data = num.linalg.norm(ddisloc, axis=1) / dt
@@ -1558,16 +1682,29 @@ def rupture_movie(
         fn_temp_path = op.join(temp_path, 'f%09d.png')
 
         for it, (t, ft) in enumerate(zip(times, fns_temp)):
+            plot_data = data[:, it]
+
+            if original_patch_indices is not None:
+                plot_data = num.zeros(original_patch_indices.shape)
+                plot_data[original_patch_indices] = data[:, it]
+
             plt = plt_base(source=source, **kwargs_base)
-            plt.draw_dynamic_data(data[:, it],
+            plt.draw_dynamic_data(plot_data,
                                   **kwargs_plt)
             plt.draw_time_contour(store, clevel=[t])
 
             plt.save(ft)
 
-        render_movie(fn_temp_path,
-                     output_path=op.join(
-                         fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
+        fn_mp4 = op.join(temp_path, 'movie.mp4')
+        render_movie(fn_temp_path, output_path=fn_mp4)
+
+        if render_as_gif:
+            render_gif(fn=fn_mp4, output_path=op.join(
+                fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
+                loops=gif_loops)
+        else:
+            shutil.move(fn_mp4, op.join(
+                fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
 
         if store_images:
             fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))

src/scenario/sources/pseudodynrupture.py

@@ -7,7 +7,7 @@ from __future__ import absolute_import, division, print_function
 import numpy as num
 
 from pyrocko.guts import Float, Int
-from pyrocko import moment_tensor, gf
+from pyrocko import gf
 
 from .base import SourceGenerator
 

test/base/test_ahfull.py

@@ -142,7 +142,7 @@ class AhfullTestCase(unittest.TestCase):
                 'displacement',
                 s.deltat, 0.,
                 out_x, out_y, out_z,
-                ahfullgreen.Gauss(s.tau))
+                ahfullgreen.AhfullgreenSTFGauss(tau=s.tau))
 
             trs = []
             for out, comp in zip([out_x, out_y, out_z], 'NED'):

test/gf/test_gf.py

@@ -90,7 +90,7 @@ def _make_traces_homogeneous(
             material.vp, material.vs, material.rho, material.qp, material.qs,
             x, f, m6, 'displacement',
             deltat, tmin - (times[ielement] - toff), outx, outy, outz,
-            stf=ahfullgreen.Impulse())
+            stf=ahfullgreen.AhfullgreenSTFImpulse())
 
     for comp, o in zip(comps, (outx, outy, outz)):
         tr = trace.Trace(
@@ -1377,7 +1377,10 @@ class GFTestCase(unittest.TestCase):
         store_dir = mkdtemp(prefix=store_id)
         self.tempdirs.append(store_dir)
 
-        gf.store.Store.create_editables(store_dir, config=config)
+        conf = fomosto_ahfullgreen.AhfullgreenConfig()
+
+        gf.store.Store.create_editables(
+            store_dir, config=config, extra={'ahfullgreen': conf})
 
         store = gf.store.Store(store_dir, 'r')
         store.make_ttt()

test/gf/test_gf_qseis.py

@@ -239,6 +239,8 @@ mantle
         store_id_qseis = 'homogeneous_qseis'
         store_id_ahfull = 'homogeneous_ahfull'
 
+        ahconf = ahfullgreen.AhfullgreenConfig()
+
         qsconf = qseis.QSeisConfig()
         qsconf.qseis_version = '2006a'
 
@@ -338,7 +340,7 @@ mantle
         self.tempdirs.append(store_dir_ahfull)
 
         gf.store.Store.create_editables(
-            store_dir_ahfull, config=config)
+            store_dir_ahfull, config=config, extra={'ahfullgreen': ahconf})
 
         store = gf.store.Store(store_dir_ahfull, 'r')
         store.make_ttt()

test/gf/test_gf_source_types.py

@@ -82,9 +82,9 @@ class GFSourceTypesTestCase(unittest.TestCase):
             nucleation_y=0.)
 
         pdr.tractions = gf.tractions.HomogeneousTractions(
-            t_strike=0.,
-            t_dip=0.,
-            t_normal=-.5e6)
+            strike=0.,
+            dip=0.,
+            normal=-.5e6)
 
         time = num.max(pdr.get_patch_attribute('time')) * 0.5
 

test/gf/test_gf_sources.py

@@ -269,6 +269,9 @@ class GFSourcesTestCase(unittest.TestCase):
             if not hasattr(S, 'discretize_basesource'):
                 continue
 
+            if S.T.classname == 'PseudoDynamicRupture':
+                continue
+
             for t in [0.0, util.str_to_time('2014-01-01 10:00:00')]:
                 source = default_source(S, time=t)
                 dsource = source.discretize_basesource(

test/gf/test_gf_static.py

@@ -93,7 +93,7 @@ mantle
             source_depth_max=20.*km,
             source_depth_delta=0.5*km,
             distance_min=0.*km,
-            distance_max=40.*km,
+            distance_max=70.*km,
             distance_delta=0.5*km,
             modelling_code_id='psgrn_pscmp.%s' % version,
             earthmodel_1d=mod,
@@ -287,9 +287,9 @@ mantle
             dyn_rupture = gf.PseudoDynamicRupture(
                 nx=nx, ny=ny,
                 tractions=gf.tractions.HomogeneousTractions(
-                    t_strike=1.e4,
-                    t_dip=0.4e4,
-                    t_normal=0.1),
+                    strike=1.e4,
+                    dip=0.4e4,
+                    normal=0.1),
                 north_shift=2*km,
                 east_shift=2*km,
                 depth=6.5*km,
@@ -310,7 +310,8 @@ mantle
             t = time.time()
             # dyn_rupture.discretize_patches(store)
             engine.process(dyn_rupture, static_target)
-            map = RuptureMap(source=dyn_rupture, lat=0., lon=0., radius=40*km)
+            map = RuptureMap(source=dyn_rupture, lat=0., lon=0., radius=40*km,
+                             width=20., height=20.)
             map.draw_patch_parameter('traction')
             map.save('/tmp/test.pdf')
             return dyn_rupture.nx*dyn_rupture.ny, time.time() - t

test/gf/test_modelling_okada.py

@@ -166,6 +166,7 @@ class OkadaTestCase(unittest.TestCase):
             num.linalg.inv(num.dot(res1.T, res1))
 
         res = res1.copy()
+
         @benchmark.labeled('mat_inversion_sparse')
         def inv_sparse():
             res[res < 1e-2] = 0.