"""
Module for interseismic models.
Block-backslip model
--------------------
The fault is assumed to be locked above a certain depth "locking_depth" and
it is creeping with the rate of the defined plate- which is handled as a
rigid block.
STILL EXPERIMENTAL!
References
==========
Savage & Prescott 1978
Metzger et al. 2011
"""
from beat import utility
from beat.heart import geo_synthetics
import numpy as num
import logging
import copy
from pyrocko.orthodrome import latlon_to_ne_numpy, latlon_to_xyz, earthradius
from pyrocko.gf import RectangularSource as RS
from matplotlib.path import Path
logger = logging.getLogger('interseismic')
km = 1000.
d2r = num.pi / 180.
r2d = 180. / num.pi
non_source = set(['amplitude', 'azimuth', 'locking_depth'])
__all__ = ['geo_backslip_synthetics']
def block_mask(easts, norths, sources, east_ref, north_ref):
"""
Determine stable and moving observation points dependend on the input
fault orientation.
Parameters
----------
easts : :class:`numpy.ndarray`
east - local coordinates [m] of observations
norths : :class:`numpy.ndarray`
north - local coordinates [m] of observations
sources : list
of :class:`RectangularSource`
east_ref : float
east local coordinate [m] of stable reference
north_ref : float
north local coordinate [m] of stable reference
Returns
-------
:class:`numpy.ndarray` with zeros at stable points, ones at moving points
"""
def get_vertex(outlines, i, j):
f1 = outlines[i]
f2 = outlines[j]
print(f1, f2)
return utility.line_intersect(f1[0, :], f1[1, :], f2[0, :], f2[1, :])
tol = 2. * km
Eline = RS(
east_shift=easts.max() + tol, north_shift=0.,
strike=0., dip=90., length=1 * km)
Nline = RS(
east_shift=0., north_shift=norths.max() + tol,
strike=90, dip=90., length=1 * km)
Sline = RS(
east_shift=0., north_shift=norths.min() - tol,
strike=90, dip=90., length=1 * km)
frame = [Nline, Eline, Sline]
# collect frame lines
outlines = []
for source in sources + frame:
outline = source.outline(cs='xy')
outlines.append(utility.swap_columns(outline, 0, 1)[0:2, :])
# get polygon vertices
poly_vertices = []
for i in range(len(outlines) - 1):
poly_vertices.append(get_vertex(outlines, i, i + 1))
else:
poly_vertices.append(get_vertex(outlines, 0, -1))
print(poly_vertices, outlines)
polygon = Path(num.vstack(poly_vertices), closed=True)
ens = num.vstack([easts.flatten(), norths.flatten()]).T
ref_en = num.array([east_ref, north_ref]).flatten()
print(ens)
mask = polygon.contains_points(ens)
if not polygon.contains_point(ref_en):
return mask
else:
return num.logical_not(mask)
def block_geometry(lons, lats, sources, reference):
"""
Construct block geometry determine stable and moving parts dependend
on the reference location.
Parameters
----------
lons : :class:`num.ndarray`
Longitudes [deg] of observation points
lats : :class:`num.ndarray`
Latitudes [deg] of observation points
sources : list
of RectangularFault objects
reference : :class:`heart.ReferenceLocation`
reference location that determines the stable block
Returns
-------
:class:`num.ndarray`
mask with zeros/ones for stable/moving observation points, respectively
"""
norths, easts = latlon_to_ne_numpy(
reference.lat, reference.lon, lats, lons)
return block_mask(easts, norths, sources, east_ref=0., north_ref=0.)
def block_movement(bmask, amplitude, azimuth):
"""
Get block movements. Assumes one side of the model stable, therefore
the moving side is moving 2 times the given amplitude.
Parameters
----------
bmask : :class:`numpy.ndarray`
masked block determining stable and moving observation points
amplitude : float
slip [m] of the moving block
azimuth : float
azimuth-angle[deg] ergo direction of moving block towards North
Returns
-------
:class:`numpy.ndarray`
(n x 3) [North, East, Down] displacements [m]
"""
tmp = num.repeat(
bmask * 2. * float(amplitude), 3).reshape((bmask.shape[0], 3))
sv = utility.strike_vector(float(azimuth), order='NEZ')
return tmp * sv
def geo_block_synthetics(lons, lats, sources, amplitude, azimuth, reference):
"""
Block model: forward model for synthetic displacements(n,e,d) [m] caused by
a rigid moving block defined by the bounding geometry of rectangular
faults. The reference location determines the stable regions.
The amplitude and azimuth determines the amount and direction of the
moving block.
Parameters
----------
lons : :class:`num.ndarray`
Longitudes [deg] of observation points
lats : :class:`num.ndarray`
Latitudes [deg] of observation points
sources : list
of RectangularFault objects
amplitude : float
slip [m] of the moving block
azimuth : float
azimuth-angle[deg] ergo direction of moving block towards North
reference : :class:`heart.ReferenceLocation`
reference location that determines the stable block
Returns
-------
:class:`numpy.ndarray`
(n x 3) [North, East, Down] displacements [m]
"""
bmask = block_geometry(lons, lats, sources, reference)
return block_movement(bmask, amplitude, azimuth)
def backslip_params(azimuth, strike, dip, amplitude, locking_depth):
"""
Transforms the interseismic blockmodel parameters to fault input parameters
for the backslip model.
Parameters
----------
azimuth : float
azimuth [deg] of the block-motion towards the North
strike : float
strike-angle[deg] of the backslipping fault
dip : float
dip-angle[deg] of the back-slipping fault
amplitude : float
slip rate of the blockmodel [m/yr]
locking_depth : float
locking depth [km] of the fault
Returns
-------
dict of parameters for the back-slipping RectangularSource
"""
if dip == 0.:
raise ValueError('Dip must not be zero!')
az_vec = utility.strike_vector(azimuth)
strike_vec = utility.strike_vector(strike)
alpha = num.arccos(az_vec.dot(strike_vec))
alphad = alpha * r2d
sdip = num.sin(dip * d2r)
# assuming dip-slip is zero --> strike slip = slip
slip = num.abs(amplitude * num.cos(alpha))
opening = -amplitude * num.sin(alpha) * sdip
if alphad < 90. and alphad >= 0.:
rake = 0.
elif alphad >= 90. and alphad <= 180.:
rake = 180.
else:
raise Exception('Angle between vectors inconsistent!')
width = locking_depth * km / sdip
return dict(
slip=float(slip), opening=float(opening), width=float(width),
depth=0., rake=float(rake))
[docs]def geo_backslip_synthetics(
engine, sources, targets, lons, lats, reference,
amplitude, azimuth, locking_depth):
"""
Interseismic backslip model: forward model for synthetic
displacements(n,e,d) [m] caused by a rigid moving block defined by the
bounding geometry of rectangular faults. The reference location determines
the stable regions. The amplitude and azimuth determines the amount and
direction of the moving block.
Based on this block-movement the upper part of the crust that is not locked
is assumed to slip back. Thus the final synthetics are the superposition
of the block-movement and the backslip.
Parameters
----------
engine : :class:`pyrocko.gf.seismosizer.LocalEngine`
sources : list
of :class:`pyrocko.gf.seismosizer.RectangularSource`
Sources to calculate synthetics for
targets : list
of :class:`pyrocko.gf.targets.StaticTarget`
lons : list of floats, or :class:`numpy.ndarray`
longitudes [deg] of observation points
lats : list of floats, or :class:`numpy.ndarray`
latitudes [deg] of observation points
amplitude : float
slip [m] of the moving block
azimuth : float
azimuth-angle[deg] ergo direction of moving block towards North
locking_depth : :class:`numpy.ndarray`
locking_depth [km] of the fault(s) below there is no movement
reference : :class:`heart.ReferenceLocation`
reference location that determines the stable block
Returns
-------
:class:`numpy.ndarray`
(n x 3) [North, East, Down] displacements [m]
"""
disp_block = geo_block_synthetics(
lons, lats, sources, amplitude, azimuth, reference)
for source, ld in zip(sources, locking_depth):
source_params = backslip_params(
azimuth=azimuth, amplitude=amplitude, locking_depth=ld,
strike=source.strike, dip=source.dip)
source.update(**source_params)
disp_block += geo_synthetics(
engine=engine, targets=targets, sources=sources,
outmode='stacked_array')
return disp_block
def seperate_point(point):
"""
Seperate point into source object related components and the rest.
"""
tpoint = copy.deepcopy(point)
interseismic_point = {}
for var in non_source:
if var in tpoint.keys():
interseismic_point[var] = tpoint.pop(var)
return tpoint, interseismic_point
