Forward modeling synthetic seismograms and displacements

Calculate synthetic seismograms from a local GF store

It is assumed that a Store with store ID crust2_dd has been downloaded in advance. A list of currently available stores can be found at https://greens-mill.pyrocko.org as well as how to download such stores.

Further API documentation for the utilized objects can be found at Target, LocalEngine and DCSource.

Download gf_forward_example1.py

import os

from pyrocko.gf import LocalEngine, Target, DCSource, ws
from pyrocko import trace
from pyrocko.marker import PhaseMarker

# The store we are going extract data from:
store_id = 'iceland_reg_v2'

# First, download a Greens Functions store. If you already have one that you
# would like to use, you can skip this step and point the *store_superdirs* in
# the next step to that directory.

if not os.path.exists(store_id):
    ws.download_gf_store(site='kinherd', store_id=store_id)

# We need a pyrocko.gf.Engine object which provides us with the traces
# extracted from the store. In this case we are going to use a local
# engine since we are going to query a local store.
engine = LocalEngine(store_superdirs=['.'])

# Define a list of pyrocko.gf.Target objects, representing the recording
# devices. In this case one station with a three component sensor will
# serve fine for demonstation.
channel_codes = 'ENZ'
targets = [
    Target(
        lat=10.,
        lon=10.,
        store_id=store_id,
        codes=('', 'STA', '', channel_code))
    for channel_code in channel_codes]

# Let's use a double couple source representation.
source_dc = DCSource(
    lat=11.,
    lon=11.,
    depth=10000.,
    strike=20.,
    dip=40.,
    rake=60.,
    magnitude=4.)

# Processing that data will return a pyrocko.gf.Reponse object.
response = engine.process(source_dc, targets)

# This will return a list of the requested traces:
synthetic_traces = response.pyrocko_traces()

# In addition to that it is also possible to extract interpolated travel times
# of phases which have been defined in the store's config file.
store = engine.get_store(store_id)

markers = []
for t in targets:
    dist = t.distance_to(source_dc)
    depth = source_dc.depth
    arrival_time = store.t('begin', (depth, dist))
    m = PhaseMarker(tmin=arrival_time,
                    tmax=arrival_time,
                    phasename='P',
                    nslc_ids=(t.codes,))
    markers.append(m)

# Finally, let's scrutinize these traces.
trace.snuffle(synthetic_traces, markers=markers)
Synthetic seismograms calculated through pyrocko.gf

Synthetic seismograms calculated through pyrocko.gf displayed in Snuffler - seismogram browser and workbench. The three traces show the east, north and vertical synthetical displacement stimulated by a double-couple source at 155 km distance.

Calculate spatial surface displacement from a local GF store

In this example we create a RectangularSource and compute the spatial static/geodetic displacement caused by that rupture.

We will utilize LocalEngine, StaticTarget and SatelliteTarget.

Static displacement from a strike-slip fault calculated through pyrocko

Synthetic surface displacement from a vertical strike-slip fault, with a N104W azimuth, in the Line-of-sight (LOS), east, north and vertical directions. LOS as for Envisat satellite (Look Angle: 23., Heading:-76). Positive motion toward the satellite.

Download gf_forward_example2.py

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

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

# Setup the LocalEngine and point it to the fomosto store you just downloaded.
# *store_superdirs* is a list of directories where to look for GF Stores.
engine = gf.LocalEngine(store_superdirs=['.'])

# We define an extended source, in this case a rectangular geometry
# Centroid UTM position is defined relatively to geographical lat, lon position
# Purely lef-lateral strike-slip fault with an N104W azimuth.

km = 1e3  # for convenience

rect_source = gf.RectangularSource(
    lat=0., lon=0.,
    north_shift=0., east_shift=0., depth=6.5*km,
    width=5*km, length=8*km,
    dip=90., rake=0., strike=104.,
    slip=1.)

# We will define a grid of targets
# number in east and north directions, and total
ngrid = 80

# extension from origin in all directions
obs_size = 20.*km
ntargets = ngrid**2

# make regular line vector
norths = num.linspace(-obs_size, obs_size, ngrid)
easts = num.linspace(-obs_size, obs_size, ngrid)

# make regular grid
norths2d = num.repeat(norths, len(easts))
easts2d = num.tile(easts, len(norths))

# We initialize the satellite target and set the line of sight vectors
# direction, example of the Envisat satellite
look = 23.     # angle between the LOS and the vertical
heading = -76  # angle between the azimuth and the east (anti-clock)
theta = num.empty(ntargets)  # vertical LOS from horizontal
theta.fill(num.deg2rad(90. - look))
phi = num.empty(ntargets)  # horizontal LOS from E in anti-clokwise rotation
phi.fill(num.deg2rad(-90-heading))

satellite_target = gf.SatelliteTarget(
    north_shifts=norths2d,
    east_shifts=easts2d,
    tsnapshot=24. * 3600.,  # one day
    interpolation='nearest_neighbor',
    phi=phi,
    theta=theta,
    store_id=store_id)

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


def plot_static_los_result(result, target=0):
    '''Helper function for plotting the displacement'''

    import matplotlib.pyplot as plt

    # get target coordinates and displacements from results
    N = result.request.targets[target].coords5[:, 2]
    E = result.request.targets[target].coords5[:, 3]
    synth_disp = result.results_list[0][target].result

    # get the component names of displacements
    components = synth_disp.keys()
    fig, _ = plt.subplots(int(len(components)/2), int(len(components)/2))

    vranges = [(synth_disp[k].max(),
                synth_disp[k].min()) for k in components]

    for comp, ax, vrange in zip(components, fig.axes, vranges):

        lmax = num.abs([num.min(vrange), num.max(vrange)]).max()

        # plot displacements at targets as colored points
        cmap = ax.scatter(E, N, c=synth_disp[comp], s=10., marker='s',
                          edgecolor='face', cmap='seismic',
                          vmin=-1.5*lmax, vmax=1.5*lmax)

        ax.set_title(comp+' [m]')
        ax.set_aspect('equal')
        ax.set_xlim(-obs_size, obs_size)
        ax.set_ylim(-obs_size, obs_size)
        # We plot the modeled fault
        n, e = rect_source.outline(cs='xy').T
        ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)

        fig.colorbar(cmap, ax=ax, aspect=5)

    plt.show()


plot_static_los_result(result)

Calculate forward model of thrust event and display wrapped phase

In this example we compare the synthetic unwappred and wrapped LOS displacements caused by a thrust rupture.

Static displacement from a thrust fault calculated through pyrocko

Synthetic LOS displacements from a south-dipping thrust fault. LOS as for Sentinel-1 satellite (Look Angle: 36., Heading:-76). Positive motion toward the satellite. Left: unwrapped phase. Right: Wrapped phase.

Download gf_forward_example3.py

import os.path
import numpy as num
import matplotlib.pyplot as plt
from pyrocko import gf

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

# Ignite the LocalEngine and point it to your fomosto store, e.g. stored on a
# USB stick, which for example has the id 'Abruzzo_Ameri_static_nearfield'
engine = gf.LocalEngine(store_superdirs=['.'])

# We want to reproduce the USGS Solution of an event, e.g.
dep, strike, dip, leng, wid, rake, slip = 10.5, 90., 40., 10., 10., 90., .5

km = 1e3    # distance in kilometer, for convenienve

# We compute the magnitude of the event
potency = leng*km*wid*km*slip
rigidity = 31.5e9
m0 = potency*rigidity
mw = (2./3) * num.log10(m0) - 6.07

# We define an extended rectangular source
thrust = gf.RectangularSource(
    north_shift=0., east_shift=0.,
    depth=dep*km, width=wid*km, length=leng*km,
    dip=dip, rake=rake, strike=strike,
    slip=slip)

# We will define a grid of targets
# number in east and north directions, and total
ngrid = 40
# ngrid = 90  # for better resolution

# extension from origin in all directions
obs_size = 20.*km
ntargets = ngrid**2

# make regular line vector
norths = num.linspace(-obs_size, obs_size, ngrid)
easts = num.linspace(-obs_size, obs_size, ngrid)

# make regular grid
norths2d = num.repeat(norths, len(easts))
easts2d = num.tile(easts, len(norths))

# 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=norths2d,
    east_shifts=easts2d,
    tsnapshot=24.*3600.,    # one day
    interpolation='nearest_neighbor',
    phi=phi,
    theta=theta,
    store_id=store_id)

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

# get synthetic displacements and target coordinates from engine's 'result'
# of the first target (i target=0)
i_target = 0
N = result.request.targets[i_target].coords5[:, 2]
E = result.request.targets[i_target].coords5[:, 3]
synth_disp = result.results_list[0][i_target].result

# we get the fault projection to the surface for plotting
n, e = thrust.outline(cs='xy').T

fig, _ = plt.subplots(1, 2, figsize=(8, 4))
fig.suptitle(
    "fault: dep={:0.2f}, l={}, w={:0.2f},str={},"
    "rake={}, dip={}, slip={}, Mw={:0.3f}\n"
    "satellite: heading={}, look angle={}"
    .format(dep, leng, wid, strike, rake, dip, slip, heading, look, mw),
    fontsize=14,
    fontweight='bold')

# We shift the relative LOS displacements
los = synth_disp['displacement.los']
losrange = [(los.max(), los.min())]
losmax = num.abs([num.min(losrange), num.max(losrange)]).max()

# Plot unwrapped LOS displacements
ax = fig.axes[0]
cmap = ax.scatter(
    E, N, c=los,
    s=10., marker='s',
    edgecolor='face',
    cmap=plt.get_cmap('seismic'),
    vmin=-1.*losmax, vmax=1.*losmax)

ax.set_title('line-of-sight displacement [m]')
ax.set_aspect('equal')
ax.set_xlim(-obs_size, obs_size)
ax.set_ylim(-obs_size, obs_size)
# plot fault outline
ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
# We underline the tip of the thrust
ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)

fig.colorbar(cmap, ax=ax, orientation='vertical', aspect=5, shrink=0.5)

# We simulate a C-band interferogram for this source
c_lambda = 0.056
insar_phase = -num.mod(los, c_lambda/2.)/(c_lambda/2.)*2.*num.pi - num.pi

# We plot wrapped phase
ax = fig.axes[1]
cmap = ax.scatter(
    E, N, c=insar_phase,
    s=10., marker='s',
    edgecolor='face',
    cmap=plt.get_cmap('gist_rainbow'))

ax.set_xlim(-obs_size, obs_size)
ax.set_ylim(-obs_size, obs_size)
ax.set_title('simulated interferogram')
ax.set_aspect('equal')

# plot fault outline
ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
# We outline the top edge of the fault with a thick line
ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)

fig.colorbar(cmap, orientation='vertical', shrink=0.5, aspect=5)
plt.show()

Combining severals sources

In this example we combine two rectangular sources and plot the forward model in profile.

../../_images/gf_static_several.png

Synthetic LOS displacements from a flower-structure made of one strike-slip fault and one thrust fault. LOS as for Sentinel-1 satellite (Look Angle: 36., Heading:-76). Positive motion toward the satellite.

Download 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=num.float)
            lons = num.array(
                [subsource.lon for subsource in subsources], dtype=num.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):
    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.)