gf_forward_example4.py

import os.path
import numpy as num
from pyrocko import gf
from pyrocko.guts import List


class CombiSource(gf.Source):
    '''Composite source model.'''

    discretized_source_class = gf.DiscretizedMTSource

    subsources = List.T(gf.Source.T())

    def __init__(self, subsources=[], **kwargs):

        if subsources:

            lats = num.array(
                [subsource.lat for subsource in subsources], dtype=float)
            lons = num.array(
                [subsource.lon for subsource in subsources], dtype=float)

            assert num.all(lats == lats[0]) and num.all(lons == lons[0])
            lat, lon = lats[0], lons[0]

            # if not same use:
            # lat, lon = center_latlon(subsources)

            depth = float(num.mean([p.depth for p in subsources]))
            t = float(num.mean([p.time for p in subsources]))
            kwargs.update(time=t, lat=float(lat), lon=float(lon), depth=depth)

        gf.Source.__init__(self, subsources=subsources, **kwargs)

    def get_factor(self):
        return 1.0

    def discretize_basesource(self, store, target=None):
        dsources = []
        t0 = self.subsources[0].time
        for sf in self.subsources:
            assert t0 == sf.time
            ds = sf.discretize_basesource(store, target)
            ds.m6s *= sf.get_factor()
            dsources.append(ds)

        return gf.DiscretizedMTSource.combine(dsources)


# Download a Greens Functions store, programmatically.
store_id = 'gf_abruzzo_nearfield_vmod_Ameri'
if not os.path.exists(store_id):
    gf.ws.download_gf_store(site='kinherd', store_id=store_id)

km = 1e3   # distance in kilometer

# We define a grid for the targets.
left, right, bottom, top = -10*km, 10*km, -10*km, 10*km
ntargets = 1000

# Ignite the LocalEngine and point it to fomosto stores stored on a
# USB stick, for this example we use a static store with id 'static_store'
engine = gf.LocalEngine(store_superdirs=['.'])
store_id = 'gf_abruzzo_nearfield_vmod_Ameri'

# We define two finite sources
# The first one is a purely vertical strike-slip fault
strikeslip = gf.RectangularSource(
    north_shift=0., east_shift=0.,
    depth=6*km, width=4*km, length=10*km,
    dip=90., rake=0., strike=90.,
    slip=1.)

# The second one is a ramp connecting to the root of the strike-slip fault
# ramp north shift (n) and width (w) depend on its dip angle and on
# the strike slip fault width
n, w = 2/num.tan(num.deg2rad(45.)), 2.*(2./(num.sin(num.deg2rad(45.))))
thrust = gf.RectangularSource(
    north_shift=n*km, east_shift=0.,
    depth=6*km, width=w*km, length=10*km,
    dip=45., rake=90., strike=90.,
    slip=0.5)

# We initialize the satellite target and set the line of site vectors
# Case example of the Sentinel-1 satellite:
# Heading: -166 (anti clockwise rotation from east)
# Average Look Angle: 36 (from vertical)
heading = -76
look = 36.
phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clockwise rotation
theta = num.empty(ntargets)  # Vertical LOS from horizontal
phi.fill(num.deg2rad(-90. - heading))
theta.fill(num.deg2rad(90. - look))

satellite_target = gf.SatelliteTarget(
    north_shifts=num.random.uniform(bottom, top, ntargets),
    east_shifts=num.random.uniform(left, right, ntargets),
    tsnapshot=24.*3600.,
    interpolation='nearest_neighbor',
    phi=phi,
    theta=theta,
    store_id=store_id)

# We combine the two sources here
patches = [strikeslip, thrust]
sources = CombiSource(subsources=patches)

# The computation is performed by calling process on the engine
result = engine.process(sources, [satellite_target])


def plot_static_los_profile(result, strike, l, w, x0, y0):  # noqa
    import matplotlib.pyplot as plt
    import matplotlib.cm as cm
    import matplotlib.colors as mcolors
    fig, _ = plt.subplots(1, 2, figsize=(8, 4))

    # strike,l,w,x0,y0: strike, length, width, x, and y position
    # of the profile
    strike = num.deg2rad(strike)
    # We define the parallel and perpendicular vectors to the profile
    s = [num.sin(strike), num.cos(strike)]
    n = [num.cos(strike), -num.sin(strike)]

    # We define the boundaries of the profile
    ypmax, ypmin = l/2, -l/2
    xpmax, xpmin = w/2, -w/2

    # We define the corners of the profile
    xpro, ypro = num.zeros((7)), num.zeros((7))
    xpro[:] = x0-w/2*s[0]-l/2*n[0], x0+w/2*s[0]-l/2*n[0], \
        x0+w/2*s[0]+l/2*n[0], x0-w/2*s[0]+l/2*n[0], x0-w/2*s[0]-l/2*n[0], \
        x0-l/2*n[0], x0+l/2*n[0]

    ypro[:] = y0-w/2*s[1]-l/2*n[1], y0+w/2*s[1]-l/2*n[1], \
        y0+w/2*s[1]+l/2*n[1], y0-w/2*s[1]+l/2*n[1], y0-w/2*s[1]-l/2*n[1], \
        y0-l/2*n[1], y0+l/2*n[1]

    # We get the forward model from the engine
    N = result.request.targets[0].coords5[:, 2]
    E = result.request.targets[0].coords5[:, 3]
    result = result.results_list[0][0].result

    # We first plot the surface displacements in map view
    ax = fig.axes[0]
    los = result['displacement.los']
    levels = num.linspace(los.min(), los.max(), 50)

    cmap = ax.tricontourf(E, N, los, cmap=plt.get_cmap('seismic'),
                          levels=levels)

    for sourcess in patches:
        fn, fe = sourcess.outline(cs='xy').T
        ax.fill(fe, fn, color=(0.5, 0.5, 0.5), alpha=0.5)
        ax.plot(fe[:2], fn[:2], linewidth=2., color='black', alpha=0.5)

    # We plot the limits of the profile in map view
    ax.plot(xpro[:], ypro[:], color='black', lw=1.)
    # plot colorbar
    fig.colorbar(cmap, ax=ax, orientation='vertical', aspect=5)
    ax.set_title('Map view')
    ax.set_aspect('equal')

    # We plot displacements in profile
    ax = fig.axes[1]
    # We compute the perpendicular and parallel components in the profile basis
    yp = (E-x0)*n[0]+(N-y0)*n[1]
    xp = (E-x0)*s[0]+(N-y0)*s[1]
    los = result['displacement.los']

    # We select data encompassing the profile
    index = num.nonzero(
        (xp > xpmax) | (xp < xpmin) | (yp > ypmax) | (yp < ypmin))

    ypp, losp = num.delete(yp, index), \
        num.delete(los, index)

    # We associate the same color scale to the scatter plot
    norm = mcolors.Normalize(vmin=los.min(), vmax=los.max())
    m = cm.ScalarMappable(norm=norm, cmap=plt.get_cmap('seismic'))
    facelos = m.to_rgba(losp)
    ax.scatter(
        ypp, losp,
        s=0.3, marker='o', color=facelos, label='LOS displacements')

    ax.legend(loc='best')
    ax.set_title('Profile')

    plt.show()


plot_static_los_profile(result, 110., 18*km, 5*km, 0., 0.)