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

- bugfixes of buggy magnitude scaling (fix with new rescale_slip method)
- bugfixes plotting
- changed handling of Lamee parameters for better future implementation
	possibilities of heterogeneous linear coefficient calculators
- Update gmtpy import
feature/orthodrom-crossdist
mmetz 2 years ago committed by Sebastian Heimann
parent 3c2c27bf57
commit 5ae55b9155

@ -23,6 +23,25 @@ Download :download:`gf_forward_example1.py </../../examples/gf_forward_example1.
Synthetic seismograms calculated through :class:`pyrocko.gf` displayed in :doc:`/apps/snuffler/index`. The three traces show the east, north and vertical synthetical displacement stimulated by a double-couple source at 155 km distance.
Calculate synthetic seismograms using the Pseudo Dynamic Rupture
----------------------------------------------------------------
Download :download:`gf_forward_pseudo_rupture_waveforms.py </../../examples/gf_forward_pseudo_rupture_waveforms.py>`
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_waveforms.py
:language: python
.. figure :: /static/gf_forward_pseudo_rupture_waveforms.png
:align: center
:width: 90%
:alt: Synthetic seismogram calculated through pyrocko.gf using :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
Synthetic seismogram calculated through :class:`pyrocko.gf` using the
:py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
Calculate spatial surface displacement from a local GF store
------------------------------------------------------------
Shear dislocation - Earthquake
@ -76,11 +95,44 @@ Download :download:`okada_forward_example.py </../../examples/okada_forward_exam
.. literalinclude :: /../../examples/okada_forward_example.py
:language: python
.. figure :: /static/okada_forward_example.png
:align: center
:width: 90%
:alt: Surface displacements derived from a set of :py:class:`~pyrocko.modelling.okada.OkadaSource`
Surface displacements (3 components and absolute value) calculated using a
set of :py:class:`~pyrocko.modelling.okada.OkadaSource`.
.. rubric:: Footnotes
.. [#f1] Okada, Y., Gravity and potential changes due to shear and tensile faults in a half-space. In: Journal of Geophysical Research 82.2, 10181040. doi:10.1029/92JB00178, 1992.
Calculate spatial surface displacement using the Pseudo Dynamic Rupture
-----------------------------------------------------------------------
In this example we create a :class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture` and compute the spatial static displacement at the surface invoked by that rupture [#f2]_.
Download :download:`gf_forward_pseudo_rupture_static.py </../../examples/gf_forward_pseudo_rupture_static.py>`
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_static.py
:language: python
.. figure :: /static/gf_forward_pseudo_rupture_static.png
:align: center
:width: 90%
:alt: Surface displacements derived from a :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
Vertical surface displacements derived from a
:py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`. They are compared
to vertical static displacements calculated using the
:py:class:`~pyrocko.gf.seismosizer.RectangularSource`.
.. rubric:: Footnotes
.. [#f2] Okada, Y., Gravity and potential changes due to shear and tensile faults in a half-space. In: Journal of Geophysical Research 82.2, 10181040. doi:10.1029/92JB00178, 1992.
Calculate spatial surface displacement and export Kite scenes
-------------------------------------------------------------

@ -251,3 +251,81 @@ Download :download:`plot_directivity.py </../../examples/plot_directivity.py>`
emphasizing the directivity effects of the source (``envelope=True``).
Same source model: Mw 6.8 2020 Elazig-Sevrice earthquake.
Pseudo dynamic rupture - slip map, slip movie, source plots
-----------------------------------------------------------
The different attributes, rupture dislocations and their evolution over time
of the :py:class:`~pyrocko.gf.seismoszier.PseudoDynamicRupture` can be
inspected and illustrated in different ways from map view to small gifs. The
illustration of patch wise attributes is also possible with the built-in
module :py:mod:`~pyrocko.plot.dynamic_rupture`.
Maps of the given patch attributes or the rupture dislocation at any time can
be displayed using :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`.
Download :download:`dynamic_rupture_map.py
</../../examples/dynamic_rupture_map.py>`
.. literalinclude :: /../../examples/dynamic_rupture_map.py
:language: python
.. figure :: /static/dynamic_map_tractions.png
:align: center
:alt: Stress drop plotted as a map acting on the Pseudo Dynamic Rupture
Length of the stress drop vectors, which act on each subfault (patch) of
the :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
.. figure :: /static/dynamic_map_dislocation_3s.png
:align: center
:alt: Dislocation of the Pseudo Dynamic Rupture after 3 s of rupture
Shown is the length of the dislocation vectors of each subfault 3 s after
the rupture initiation. The rupture nucleation point is marked with the
red dot, the contour line indicates the tip of the rupture front. Arrows
show the length and direction of the slip vectors on the plane (shear
only).
On plane views of the given patch attributes or the rupture dislocation at any
time can be displayed using
:py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`. It also allows to
inspect single patch time lines as slip, slip rate or moment rate.
Download :download:`dynamic_rupture_viewer.py
</../../examples/dynamic_rupture_viewer.py>`
.. literalinclude :: /../../examples/dynamic_rupture_viewer.py
:language: python
.. figure :: /static/dynamic_view_source_traction.png
:align: center
:alt: Stress drop on the Pseudo Dynamic Rupture plotted as on plane view
Length of the stress drop vectors, which act on each subfault (patch) of
the :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
.. figure :: /static/dynamic_view_source_dislocation.png
:align: center
:alt: Dislocation of the Pseudo Dynamic Rupture after 3 s of rupture
Shown is the length of the dislocation vectors of each subfault 3 s after
the rupture initiation. The rupture nucleation point is marked with the
red dot, the contour line indicates the tip of the rupture front.
.. figure :: /static/dynamic_view_source_moment.png
:align: center
:alt: Cumulative seismic moment release of the Pseudo Dynamic Rupture
Shown is the sum of all subfault seismic moment releases of the
:py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
.. figure :: /static/dynamic_view_patch_moment.png
:align: center
:alt: Cumulative seismic moment of one subfault of the Rupture
Each patch has an individual time dependent moment release function. Here
the cumulative seismic moment over time is shown for the patch at 4th
position along strike and 4th position down dip.

Binary file not shown.

After

Width:  |  Height:  |  Size: 48 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 51 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 51 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 57 KiB

After

Width:  |  Height:  |  Size: 50 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 17 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 32 KiB

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 624 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 609 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 60 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 18 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 40 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 160 KiB

@ -53,9 +53,9 @@ Simple pseudo dynamic rupture forward model
-------------------------------------------
We create a simple forward model and calculate waveforms for one seismic station (:py:class:`~pyrocko.gf.targets.Target`) at about 14 km distance - The tractions will be aligned to force the defined ``rake``. The modeled waveform is displayed in the *snuffler* GUI.
Download :download:`gf_forward_pseudo_rupture_simple.py </../../examples/gf_forward_pseudo_rupture_simple.py>`
Download :download:`gf_forward_pseudo_rupture_basic.py </../../examples/gf_forward_pseudo_rupture_basic.py>`
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_simple.py
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_basic.py
:language: python
@ -92,12 +92,12 @@ Traction and rupture evolution
Initialize a simple dynamic rupture with uniform rake tractions and visualize the tractions and rupture propagation using the :py:mod:`pyrocko.plot.dynamic_rupture` module.
Download :download:`gf_forward_pseudo_rupture_simple_plot.py </../../examples/gf_forward_pseudo_rupture_simple_plot.py>`
Download :download:`gf_forward_pseudo_rupture_basic_plot.py </../../examples/gf_forward_pseudo_rupture_basic_plot.py>`
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_simple_plot.py
.. literalinclude :: /../../examples/gf_forward_pseudo_rupture_basic_plot.py
:language: python
.. figure :: /static/dynamic_simple_tractions.png
.. figure :: /static/dynamic_basic_tractions.png
:align: center
:width: 70%
:alt: Rupture propagation and tractions of a simple dynamic rupture source
@ -121,7 +121,7 @@ We can investigate the rupture propagation speed :math:`v_r` with :py:meth:`~pyr
plot.show_plot()
.. figure :: /static/dynamic_simple_vr.png
.. figure :: /static/dynamic_basic_vr.png
:align: center
:width: 70%
:alt: Rupture propagation and tractions of a simple dynamic rupture source
@ -145,7 +145,7 @@ Dislocations of the dynamic rupture source can be plotted with :py:meth:`~pyrock
plot.show_plot()
.. figure :: /static/dynamic_simple_dislocations.png
.. figure :: /static/dynamic_basic_dislocations.png
:align: center
:width: 70%
:alt: Rupture propagation and dislocation of a simple dynamic rupture source
@ -173,7 +173,7 @@ We can animate the rupture evolution using the :py:func:`pyrocko.plot.dynamic_ru
<center>
<video width="70%" controls>
<source src="https://pyrocko.org/media/dynamic_rupt_simple_dislocation.mp4" type="video/mp4">
<source src="https://data.pyrocko.org/media/dislocation_view_movie.mp4" type="video/mp4">
Your browser does not support the video tag.
</video>
</center>

@ -3,6 +3,7 @@ from pyrocko.example import get_example_data
from pyrocko import model
from pyrocko import moment_tensor as pmt
from pyrocko.plot import gmtpy
gmtpy.check_have_gmt()
# Download example data

@ -4,7 +4,7 @@ import numpy as num
from scipy.interpolate import RegularGridInterpolator as scrgi
from pyrocko.plot.automap import Map
from pyrocko import gmtpy
from pyrocko.plot import gmtpy
import pyrocko.orthodrome as otd
gmtpy.check_have_gmt()

@ -36,11 +36,11 @@ tracts = tractions.TractionComposition(
# Let's define the source now with its extension, orientation etc.
source = gf.PseudoDynamicRupture(
lat=-21.,
lon=32.,
lat=6.45,
lon=37.06,
length=30000.,
width=10000.,
strike=165.,
strike=215.,
dip=45.,
anchor='top',
gamma=0.6,
@ -59,11 +59,11 @@ source.discretize_patches(store)
# Define the rupture map object parameters as image center coordinates,
# radius in m and the size of the image
map_kwargs = dict(
lat=-21.,
lon=32.,
radius=15000.,
width=30.,
height=30.,
lat=6.50,
lon=37.06,
radius=20000.,
width=25.,
height=25.,
source=source,
show_topo=True)
@ -80,7 +80,7 @@ m.save('traction_map.png')
# contour at 3 s is displayed together with nucleation point.
m = dynamic_rupture.RuptureMap(**map_kwargs)
m.draw_dislocation(time=3, cmap='summer')
m.draw_dislocation_vector(time=3, S='i10.', I='x50')
m.draw_dislocation_vector(time=3, S='i15.', I='x20')
m.draw_time_contour(store, clevel=[3])
m.draw_dislocation_contour(time=3, clevel=[0.15])
m.draw_nucleation_point()

@ -60,17 +60,20 @@ source.discretize_patches(store)
# over time with a sampling interval of 1 s
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_source_dynamics(variable='moment', deltat=1, store=store)
viewer.save('dynamic_view_source_moment.png')
viewer.show_plot()
# Initialize the viewer and display the seismic moment release of one single
# boundary element (patch) over time
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_patch_dynamics(variable='moment', nx=3, ny=3, deltat=1)
viewer.save('dynamic_view_patch_moment.png')
viewer.show_plot()
# Initialize the viewer and display the traction vector length.
viewer.draw_patch_parameter('traction')
viewer.draw_nucleation_point()
viewer.save('dynamic_view_source_traction.png')
viewer.show_plot()
# Initialize the viewer and generate a more complex plot of the dislocation at
@ -81,4 +84,5 @@ viewer.draw_dislocation(time=3, cmap='summer')
viewer.draw_time_contour(store, clevel=[3])
viewer.draw_dislocation_contour(time=3, clevel=[0.15])
viewer.draw_nucleation_point()
viewer.save('dynamic_view_source_dislocation.png')
viewer.show_plot()

@ -1,125 +0,0 @@
import numpy as num
import os.path as op
import matplotlib.pyplot as plt
from pyrocko import gf
km = 1e3
d2r = 180./num.pi
# Download a Greens Functions store, programmatically.
store_id = 'gf_abruzzo_nearfield_vmod_Ameri'
store_id_dynamic = 'iceland_reg_v2'
if not op.exists(store_id):
gf.ws.download_gf_store(site='kinherd', store_id=store_id)
if not op.exists(store_id_dynamic):
gf.ws.download_gf_store(site='kinherd', store_id=store_id_dynamic)
engine = gf.LocalEngine(store_superdirs=['.'])
store = engine.get_store(store_id)
strike = 0.
dip = 90.
dep = 3*km
leng = 2*km
wid = 2*km
source_params = dict(
north_shift=2*km,
east_shift=2*km,
depth=dep,
width=wid,
length=leng,
dip=dip,
strike=strike,
magnitude=6.,
anchor='top',
decimation_factor=4)
dyn_rupture = gf.PseudoDynamicRupture(
nx=5, ny=5,
tractions=gf.tractions.HomogeneousTractions(
strike=1.e4,
dip=0.e4,
normal=0.),
**source_params)
dyn_rupture.discretize_patches(store)
slip = dyn_rupture.get_slip()
rake = num.arctan2(slip[:, 1].mean(), slip[:, 0].mean())
print('rake', float(rake*d2r))
depths = dyn_rupture.get_patch_attribute('depth')
rect_rupture = gf.RectangularSource(
rake=float(rake*d2r),
**source_params)
# Define a grid of targets
# number in east and north directions
# ngrid = 40
ngrid = 90 # for better resolution
# extension from origin in all directions
obs_size = 10.*km
ntargets = ngrid**2
norths = num.linspace(-obs_size, obs_size, ngrid)
easts = num.linspace(-obs_size, obs_size, ngrid)
norths2d = num.repeat(norths, len(easts))
easts2d = num.tile(easts, len(norths))
waveform_target = gf.Target(
lat=0., lon=0.,
east_shift=10*km, north_shift=0.*km,
store_id=store_id_dynamic)
static_target = gf.StaticTarget(
lats=num.zeros(norths2d.size), lons=num.zeros(norths2d.size),
north_shifts=norths2d.ravel(),
east_shifts=easts2d.ravel(),
interpolation='nearest_neighbor',
store_id=store_id)
result = engine.process(rect_rupture, static_target)
targets_static = result.request.targets_static
N = targets_static[0].coords5[:, 2]
E = targets_static[0].coords5[:, 3]
synth_disp_rect = result.results_list[0][0].result
result = engine.process(dyn_rupture, static_target)
targets_static = result.request.targets_static
N = targets_static[0].coords5[:, 2]
E = targets_static[0].coords5[:, 3]
synth_disp_dyn = result.results_list[0][0].result
down_rect = synth_disp_rect['displacement.d']
down_dyn = synth_disp_dyn['displacement.d']
down_diff = down_rect - down_dyn
print('rect', down_rect.max())
print('dyn', down_dyn.max())
fig, axes = plt.subplots(3, 1)
for ax, down in zip(axes, (down_rect, down_dyn, down_diff)):
cmap = ax.scatter(
E, N, c=down,
s=10., marker='s',
edgecolor='face',
cmap=plt.get_cmap('seismic'))
fig.colorbar(cmap, ax=ax, orientation='vertical', aspect=5, shrink=0.5)
plt.show()

@ -1,155 +0,0 @@
import logging
import os
import numpy as num
import matplotlib.pyplot as plt
from pyrocko import trace, util
from pyrocko.gf import PseudoDynamicRupture, LocalEngine, Target, tractions, ws
from pyrocko.plot import dynamic_rupture
from pyrocko.plot import mpl_graph_color
logger = logging.getLogger('pyrocko.examples.gf_forward_pseudo_rupture2')
util.setup_logging(levelname='info')
d2r = num.pi / 180.
km2m = 1000.
# Example of a self-similar traction field with increasing complexity used for
# the Pseudo Dynamic Rupture source model.
# 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=['.'], default_store_id=store_id)
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Length and Width are defined based from Wells and Coppersmith (1994).
mag = 7.5
length = 10**(-2.44 + 0.59 * mag) * km2m
width = 10**(-1.01 + 0.32 * mag) * km2m
nx, ny = int(num.ceil(length / 2000.)), int(num.ceil(width / 2000.))
nyq = int(num.floor(num.min((nx, ny))/2.))
logger.info('nx: %i, ny: %i' % (nx, ny))
logger.info('ranks <= %i for no aliasing' % nyq)
# Define the ranks (maximum wavenumber) and phases for the cosine functions
# in the SelfSimiliarTractions
ranks = num.array([1, 2, 5, 10])
phases = num.array([
-0.9120799, 1.40519485, -0.80165805, 1.65676832, 1.04067916, -2.58736667,
1.27630965, -2.55843096, 2.13857185, 1.01601178])
# Let's create the PseudoDynamicRupture using self similar tractions
source = PseudoDynamicRupture(
lat=0.,
lon=0.,
length=length,
width=width,
depth=10. * km2m,
strike=0.,
dip=0.,
anchor='top',
gamma=0.6,
nucleation_x=0.25,
nucleation_y=-0.5,
nx=nx,
ny=ny,
pure_shear=True,
smooth_rupture=True,
magnitude=mag,
tractions=tractions.SelfSimilarTractions(
rank=ranks[0], rake=0, phases=phases[:1]))
# The source needs to be discretized into finite faults (patches) with
# associated elastic parameters taken from the store.
source.discretize_patches(store)
# A key element of the PseudoDynamicRupture is the linkage of the tractions on
# the patches with their dislocations. The linear coefficients describing the
# link are obtained based on Okada (1992) and a boundary element method
source.calc_coef_mat()
synthetic_traces = []
channel_codes = 'ENZ'
for rank in ranks:
logger.info('Modelling for rank %i' % rank)
# Update source traction rank and phases for increasing number of summed
# cosines
source.tractions.rank = rank
source.tractions.phases = phases[:rank]
# Display absolut tractions and final absolut dislocations
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_patch_parameter('traction')
viewer.draw_time_contour(store)
viewer.show_plot()
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_dislocation()
viewer.draw_time_contour(store)
viewer.show_plot()
# 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.logger.info('Modelling synthetic waveforms')
targets = [
Target(
lat=3.,
lon=2.,
store_id=store_id,
codes=('', 'STA', '%02d' % rank, channel_code))
for channel_code in channel_codes]
# Processing that data will return a pyrocko.gf.Reponse object.
response = engine.process(source, targets)
# This will return a list of the requested traces:
synthetic_traces += response.pyrocko_traces()
# Finally, let's scrutinize these traces.
trace.snuffle(synthetic_traces)
# Plot the component-wise amplitude spectra
fig, axes = plt.subplots(3, 1)
linestyles = ['solid', 'dashed', 'dotted', 'dashdot', (0, (3, 1, 1, 1, 1, 1))]
for c, ax in zip(channel_codes, axes):
selected_traces = [tr for tr in synthetic_traces if tr.channel == c]
for i, (tr, linestyle) in enumerate(zip(selected_traces, linestyles)):
tr.ydata -= tr.ydata.mean()
freqs, amps = tr.spectrum(tfade=None)
amps = num.abs(amps)
ax.loglog(
freqs,
amps,
linestyle=linestyle,
c=mpl_graph_color(i),
label='ratio = %g' % int(tr.location))
ax.set_ylim((1e-5, 1.))
ax.set_title(c)
ax.set_ylabel('amplitude [counts]')
ax.legend(loc='best')
axes[-1].set_xlabel('frequency [Hz]')
plt.show()

@ -1,152 +0,0 @@
import logging
import os
import numpy as num
import matplotlib.pyplot as plt
from pyrocko import trace, util
from pyrocko.gf import PseudoDynamicRupture, LocalEngine, Target, tractions, ws
from pyrocko.plot import dynamic_rupture
from pyrocko.plot import mpl_graph_color
logger = logging.getLogger('pyrocko.examples.gf_forward_pseudo_rupture3')
util.setup_logging(levelname='info')
d2r = num.pi / 180.
km2m = 1000.
# Example of a fractal perturbed traction field used for the Pseudo Dynamic
# Rupture source model.
# 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=['.'], default_store_id=store_id)
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Length and Width are defined based from Wells and Coppersmith (1994).
mag = 7.5
length = 10**(-2.44 + 0.59 * mag) * km2m
width = 10**(-1.01 + 0.32 * mag) * km2m
nx = int(num.ceil(length / (2. * km2m)))
ny = int(num.ceil(width / (2. * km2m)))
logger.info('nx: %i, ny: %i' % (nx, ny))
# Define fractal traction to directed traction ratios for display of the effect
traction_ratio = num.power(num.ones(5) * 10, num.arange(-2, 3))
# Let's create the PseudoDynamicRupture using combined fractal and directed
# tractions
source = PseudoDynamicRupture(
lat=0.,
lon=0.,
length=length,
width=width,
depth=10. * km2m,
strike=0.,
dip=0.,
anchor='top',
gamma=0.6,
nucleation_x=0.25,
nucleation_y=-0.5,
nx=nx,
ny=ny,
pure_shear=True,
smooth_rupture=True,
magnitude=mag,
tractions=tractions.TractionComposition(components=[
tractions.FractalTractions(
rake=0, rseed=None, traction_max=traction_ratio[0]),
tractions.DirectedTractions(
rake=0, traction=1.)]))
# The source needs to be discretized into finite faults (patches) with
# associated elastic parameters taken from the store.
source.discretize_patches(store)
# A key element of the PseudoDynamicRupture is the linkage of the tractions on
# the patches with their dislocations. The linear coefficients describing the
# link are obtained based on Okada (1992) and a boundary element method
source.calc_coef_mat()
synthetic_traces = []
channel_codes = 'ENZ'
for it, t in enumerate(traction_ratio):
# Display absolut tractions and final absolut dislocations
logger.info('Modelling for fractal traction for traction ratio %g', t)
source.tractions.components[0].traction_max = t
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_patch_parameter('traction')
viewer.draw_time_contour(store)
viewer.show_plot()
viewer = dynamic_rupture.RuptureView(source=source)
viewer.draw_dislocation()
viewer.draw_time_contour(store)
viewer.show_plot()
# 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.logger.info('Modelling synthetic waveforms')
targets = [
Target(
lat=3.,
lon=2.,
store_id=store_id,
codes=('', 'STA', '%g' % it, channel_code))
for channel_code in channel_codes]
# Processing that data will return a pyrocko.gf.Reponse object.
response = engine.process(source, targets)
# This will return a list of the requested traces:
synthetic_traces += response.pyrocko_traces()
# Finally, let's scrutinize these traces.
trace.snuffle(synthetic_traces)
# Plot the component-wise amplitude spectra
fig, axes = plt.subplots(3, 1)
linestyles = ['solid', 'dashed', 'dotted', 'dashdot', (0, (3, 1, 1, 1, 1, 1))]
for c, ax in zip(channel_codes, axes):
selected_traces = [tr for tr in synthetic_traces if tr.channel == c]
for i, (tr, linestyle) in enumerate(zip(selected_traces, linestyles)):
tr.ydata -= tr.ydata.mean()
freqs, amps = tr.spectrum(tfade=None)
amps = num.abs(amps)
ax.loglog(
freqs,
amps,
linestyle=linestyle,
c=mpl_graph_color(i),
label='ratio = %g' % traction_ratio[int(tr.location)])
ax.set_ylim((1e-5, 1.))
ax.set_title(c)
ax.set_ylabel('amplitude [counts]')
ax.legend(loc='best')
axes[-1].set_xlabel('frequency [Hz]')
plt.show()

@ -0,0 +1,66 @@
import os.path as op
from pyrocko import gf
km = 1e3
# 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 op.exists(store_id):
gf.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 = gf.LocalEngine(store_superdirs=['.'])
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Let's define the source now with its extension, orientation etc.
dyn_rupture = gf.PseudoDynamicRupture(
lat=0.,
lon=0.,
north_shift=2.*km,
east_shift=2.*km,
depth=3.*km,
strike=43.,
dip=89.,
rake=88.,
length=15*km,
width=5*km,
nx=10,
ny=5,
# Relative nucleation between -1. and 1.
nucleation_x=-.6,
nucleation_y=.3,
slip=1.,
anchor='top',
# Threads used for modelling
nthreads=1,
# Force pure shear rupture
pure_shear=True)
# Recalculate slip, that rupture magnitude fits given magnitude
dyn_rupture.rescale_slip(magnitude=7.0, store=store)
# Create waveform target, where synthetic waveforms are calculated for
waveform_target = gf.Target(
lat=0.,
lon=0.,
east_shift=10*km,
north_shift=10.*km,
store_id=store_id)
# Get synthetic waveforms and display them in snuffler
result = engine.process(dyn_rupture, waveform_target)
result.snuffle()

@ -0,0 +1,62 @@
import os.path as op
from pyrocko import gf
from pyrocko.plot.dynamic_rupture import RuptureView
km = 1e3
# 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 op.exists(store_id):
gf.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 = gf.LocalEngine(store_superdirs=['.'])
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Let's define the source now with its extension, orientation etc.
dyn_rupture = gf.PseudoDynamicRupture(
# At lat 0. and lon 0. (default)
north_shift=2.*km,
east_shift=2.*km,
depth=3.*km,
strike=43.,
dip=89.,
rake=88.,
length=15*km,
width=5*km,
nx=10,
ny=5,
# Relative nucleation between -1. and 1.
nucleation_x=-.6,
nucleation_y=.3,
slip=1.,
anchor='top',
# Threads used for modelling
nthreads=1,
# Force pure shear rupture
pure_shear=True)
# Recalculate slip, that rupture magnitude fits given magnitude
dyn_rupture.rescale_slip(magnitude=7.0, store=store)
plot = RuptureView(dyn_rupture, figsize=(8, 4))
plot.draw_patch_parameter('traction')
plot.draw_time_contour(store)
plot.draw_nucleation_point()
plot.save('dynamic_basic_tractions.png')
# Alternatively plot on screen
# plot.show_plot()

@ -5,42 +5,59 @@ from pyrocko.plot.dynamic_rupture import RuptureView
km = 1e3
# Download a Greens Functions store
store_id = 'crust2_ib'
# 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 op.exists(store_id):
gf.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 = gf.LocalEngine(store_superdirs=['.'])
engine = gf.LocalEngine(store_superdirs=['.'], use_config=True)
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Initialize the source model
# Let's define the source now with its extension, orientation etc.
dyn_rupture = gf.PseudoDynamicRupture(
depth=3*km,
lat=0.,
lon=0.,
north_shift=2.*km,
east_shift=2.*km,
depth=3.*km,
strike=43.,
dip=89.,
length=26*km,
width=12*km,
length=15*km,
width=5*km,
nx=10,
ny=5,
# Relative nucleation between -1. and 1.
nucleation_x=-.6,
nucleation_y=.3,
# Relation between subsurface model s-wave velocity vs
# and rupture velocity vr
gamma=0.7,
magnitude=7.,
slip=1.,
anchor='top',
nthreads=5,
nx=30,
ny=20,
# Threads used for modelling
nthreads=1,
# Force pure shear rupture
pure_shear=True,
# Tractions are in [Pa]. However here we are using relative tractions,
# Resulting waveforms will be scaled to magnitude [Mw]
# or a maximumslip [m]
#
# slip=3.,
# Resulting waveforms will be scaled to the maximumslip [m]
tractions=gf.tractions.TractionComposition(
components=[
gf.tractions.DirectedTractions(rake=56., traction=1.e6),

@ -1,50 +0,0 @@
import os.path as op
from pyrocko import gf
from pyrocko.plot.dynamic_rupture import RuptureView, rupture_movie
km = 1e3
# Download a Greens Functions store
store_id = 'crust2_ib'
if not op.exists(store_id):
gf.ws.download_gf_store(site='kinherd', store_id=store_id)
engine = gf.LocalEngine(store_superdirs=['.'], use_config=True)
store = engine.get_store(store_id)
dyn_rupture = gf.PseudoDynamicRupture(
nx=30, ny=40,
north_shift=2*km, east_shift=2*km, depth=5*km,
strike=43., dip=90., rake=88.,
width=12*km, length=26*km,
nucleation_x=-.6, nucleation_y=.3,
gamma=0.7, slip=2., anchor='top', smooth_rupture=True,
nthreads=0)
dyn_rupture.discretize_patches(store)
plot = RuptureView(dyn_rupture, figsize=(8, 4))
plot.draw_patch_parameter('traction')
plot.draw_time_contour(store)
plot.draw_nucleation_point()
plot.save('dynamic_simple_tractions.png')
# Alternatively plot on screen
# plot.show_plot()
plot = RuptureView(dyn_rupture, figsize=(8, 4))
plot.draw_dislocation()
plot.draw_time_contour(store)
plot.draw_nucleation_point()
plot.save('dynamic_simple_dislocations.png')
plot = RuptureView(dyn_rupture, figsize=(8, 4))
plot.draw_patch_parameter('vr')
plot.draw_time_contour(store)
plot.draw_nucleation_point()
plot.save('dynamic_simple_vr.png')
rupture_movie(
dyn_rupture, store, 'dislocation',
plot_type='view', figsize=(8, 4))

@ -1,48 +0,0 @@
import os.path as op
from pyrocko import gf
km = 1e3
# Download a Greens Functions store
store_id = 'crust2_ib'
if not op.exists(store_id):
gf.ws.download_gf_store(site='kinherd', store_id=store_id)
engine = gf.LocalEngine(store_superdirs=['.'], use_config=True)
store = engine.get_store(store_id)
dyn_rupture = gf.PseudoDynamicRupture(
lat=0.,
lon=0.,
north_shift=0*km,
east_shift=0*km,
depth=3*km,
width=12*km,
length=29*km,
strike=43.,
dip=89.,
rake=88.,
# Number of discrete patches
nx=15,
ny=15,
# Relation between subsurface model s-wave velocity vs
# and rupture velocity vr
gamma=0.7,
magnitude=7.,
anchor='top')
waveform_target = gf.Target(
lat=0.,
lon=0.,
east_shift=10*km,
north_shift=10.*km,
store_id=store_id)
result = engine.process(dyn_rupture, waveform_target)
result.snuffle()

@ -1,48 +0,0 @@
import os.path as op
from pyrocko import gf
from pyrocko.plot.dynamic_rupture import RuptureView
km = 1e3
# Download a Greens Functions store
store_id = 'crust2_ib'
if not op.exists(store_id):
gf.ws.download_gf_store(site='kinherd', store_id=store_id)
engine = gf.LocalEngine(store_superdirs=['.'], use_config=True)
store = engine.get_store(store_id)
dyn_rupture = gf.PseudoDynamicRupture(
# At lat 0. and lon 0. (default)
north_shift=2*km,
east_shift=2*km,
depth=3*km,
strike=43.,
dip=89.,
rake=88.,
length=26*km,
width=12*km,
nx=30,
ny=20,
# Relative nucleation between -1. and 1.
nucleation_x=-.6,
nucleation_y=.3,
magnitude=7.,
anchor='top',
# Threads used for modelling
nthreads=0)
dyn_rupture.discretize_patches(store)
plot = RuptureView(dyn_rupture, figsize=(8, 4))
plot.draw_patch_parameter('traction')
plot.draw_time_contour(store)
plot.draw_nucleation_point()
plot.save('dynamic_simple_tractions.png')
# Alternatively plot on screen
# plot.show_plot()

@ -0,0 +1,129 @@
import numpy as num
import os.path as op
import matplotlib.pyplot as plt
from pyrocko import gf
km = 1e3
d2r = num.pi / 180.
r2d = 180. / num.pi
# The store we are going extract data from:
store_id = 'gf_abruzzo_nearfield_vmod_Ameri'
# 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 op.exists(store_id):
gf.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 = gf.LocalEngine(store_superdirs=['.'])
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Let's define the source now with its extension, orientation etc.
source_params = dict(
north_shift=2. * km,
east_shift=2. * km,
depth=3.5 * km,
length=6. * km,
width=3. * km,
strike=0.,
dip=0.,
rake=45.,
anchor='top',
decimation_factor=1)
dyn_rupture = gf.PseudoDynamicRupture(
nx=1,
ny=1,
pure_shear=True,
**source_params)
# Recalculate slip, that rupture magnitude fits given magnitude
magnitude = 6.0
dyn_rupture.rescale_slip(magnitude=magnitude, store=store)
# Get rake out of slip (can differ from traction rake!)
slip = dyn_rupture.get_slip()
source_params['rake'] = num.arctan2(slip[0, 1], slip[0, 0]) * r2d
# Create similar rectangular source model with rake derivded from slip
rect_rupture = gf.RectangularSource(
magnitude=magnitude,
**source_params)
# Define static target grid to extract the surface displacement
ngrid = 40
obs_size = 10. * km
ntargets = ngrid**2
norths = num.linspace(-obs_size, obs_size, ngrid) + \
source_params['north_shift']
easts = num.linspace(-obs_size, obs_size, ngrid) + \
source_params['east_shift']
norths2d = num.repeat(norths, len(easts))
easts2d = num.tile(easts, len(norths))
static_target = gf.StaticTarget(
lats=num.ones(norths2d.size) * dyn_rupture.effective_lat,
lons=num.ones(norths2d.size) * dyn_rupture.effective_lon,
north_shifts=norths2d.ravel(),
east_shifts=easts2d.ravel(),
interpolation='nearest_neighbor',
store_id=store_id)
# Get static surface displacements for rectangular and pseudo dynamic source
result = engine.process(rect_rupture, static_target)
targets_static = result.request.targets_static
synth_disp_rect = result.results_list[0][0].result
result = engine.process(dyn_rupture, static_target)
targets_static = result.request.targets_static
N = targets_static[0].coords5[:, 2]
E = targets_static[0].coords5[:, 3]
synth_disp_dyn = result.results_list[0][0].result
# Extract static vertical displacement and plot
down_rect = synth_disp_rect['displacement.d']
down_dyn = synth_disp_dyn['displacement.d']
down_diff = down_rect - down_dyn
vmin = num.min([down_rect, down_dyn, down_diff])
vmax = num.max([down_rect, down_dyn, down_diff])
fig, axes = plt.subplots(3, 1, sharex=True)
for ax, (down, label) in zip(
axes,
zip((down_rect, down_dyn, down_diff),
(r'$u_{Z, rect}$', r'$u_{Z, dyn}$', r'$\Delta u_{Z}$'))):
cntr = ax.tricontourf(
E, N, down, levels=14, cmap='RdBu_r',
vmin=vmin,
vmax=vmax)
cbar = fig.colorbar(
cntr,
ax=ax,
orientation='vertical',
aspect=10,
shrink=1.)
cbar.ax.set_ylabel(label + ' [m]')
ax.set_ylabel('Easting [m]')
axes[-1].set_xlabel('Northing [m]')
plt.show()

@ -0,0 +1,65 @@
import os.path as op
from pyrocko import gf
km = 1e3
# The store we are going extract data from:
store_id = 'crust2_mf'
# 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 op.exists(store_id):
gf.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 = gf.LocalEngine(store_superdirs=['.'], use_config=True)
# The dynamic parameter used for discretization of the PseudoDynamicRupture are
# extracted from the stores config file.
store = engine.get_store(store_id)
# Let's define the source now with its extension, orientation etc.
dyn_rupture = gf.PseudoDynamicRupture(
# At lat 0. and lon 0. (default)
north_shift=2.*km,
east_shift=2.*km,
depth=3.*km,
strike=43.,
dip=89.,
rake=88.,
length=26.*km,
width=12.*km,
nx=10,
ny=5,
# Relative nucleation between -1. and 1.
nucleation_x=-.6,
nucleation_y=.3,
slip=1.,
anchor='top',
# Threads used for modelling
nthreads=1,
# Force pure shear rupture
pure_shear=True)
# Recalculate slip, that rupture magnitude fits given magnitude
dyn_rupture.rescale_slip(magnitude=7.0, store=store)
# Model waveforms for a single station target
waveform_target = gf.Target(
lat=0.,
lon=0.,
east_shift=10.*km,
north_shift=10.*km,
store_id=store_id)
result = engine.process(dyn_rupture, waveform_target)
result.snuffle()

@ -1,7 +1,11 @@
import numpy as num
from matplotlib import pyplot as plt
from matplotlib.ticker import FuncFormatter
from pyrocko.modelling import OkadaSource, okada_ext
from pyrocko.plot import dislocation as displt
from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
d2r = num.pi / 180.
km = 1000.
@ -47,5 +51,169 @@ result = okada_ext.okada(
source.lamb, source.shearmod, nthreads=0, rotate_sdn=False,
stack_sources=True)
def draw(
axes,
dislocation,
coordinates,
xlims=[],
ylims=[],
zero_center=False,
*args,
**kwargs):
'''
Do scatterplot of dislocation array
:param axes: container for figure elements, as plot, coordinate system etc.
:type axes: :py:class:`matplotlib.axes`
:param dislocation: Dislocation array [m]
:type dislocation: :py:class:`numpy.ndarray`, ``(N,)``
:param xlims: x limits of the plot [m]
:type xlims: optional, :py:class:`numpy.ndarray`, ``