Stress inversion of the Geysers focal mechanisms from NCEDC
Cite the NCEDC: “NCEDC (2014), Northern California Earthquake Data Center. UC Berkeley Seismological Laboratory. Dataset. doi:10.7932/NCEDC.”
Acknowledge the NCEDC: “Waveform data, metadata, or data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC.”
You may have to install mplstereonet https://github.com/joferkington/mplstereonet, and colorcet https://colorcet.holoviz.org/.
[1]:
import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import string
import cartopy as ctp
import colorcet as cc
from obspy.imaging import beachball as obb
from matplotlib.colors import Normalize
from matplotlib.cm import ScalarMappable
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from shapely import geometry
import mplstereonet as mpl
import ILSI
from time import time as give_time
# plot data in Kaverina diagram
sys.path.append(
os.path.join(
"/home/eric/WORK/software/FMC/"
)
)
import plotFMC
import functionsFMC
from functionsFMC import kave, mecclass
# set plotting parameters
import seaborn as sns
sns.set_theme(font_scale=1.3)
sns.set_style("ticks")
sns.set_palette("colorblind")
plt.rcParams["savefig.dpi"] = 200
plt.rcParams["svg.fonttype"] = "none"
# define the color palette
_colors_ = ["C0", "C2", "C1", "C4"]
[2]:
# %config InlineBackend.figure_formats = ["svg"]
[3]:
# path variables
data_path = "../data"
dataset_fn = "Geysers_dataset.csv"
# program parameter(s)
save_figures = False
# load dataset
data = pd.read_csv(os.path.join(data_path, dataset_fn), sep="\t", index_col=0, header=0)
data["date"] = pd.to_datetime(data["date"], format="%Y-%m-%d")
data
[3]:
| date | hhmm | sec | event_id | mag | lat | lon | depth | strike | dip | rake | err_strike | err_dip | err_rake | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1221 | 2010-12-03 | 1049 | 44.91 | 71046544.0 | 1.03 | 38.842833 | -122.827333 | 1.68 | 10.0 | 60 | -120.0 | 35 | 33 | 50 |
| 1223 | 2010-12-04 | 439 | 59.97 | 71492300.0 | 1.69 | 38.843000 | -122.822500 | 1.94 | 5.0 | 75 | -90.0 | 43 | 28 | 40 |
| 1228 | 2010-12-05 | 132 | 44.55 | 71492590.0 | 1.32 | 38.843000 | -122.828667 | 1.84 | 10.0 | 75 | -150.0 | 25 | 35 | 40 |
| 1229 | 2010-12-05 | 1420 | 35.92 | 71492810.0 | 1.13 | 38.840500 | -122.826667 | 1.87 | 55.0 | 70 | -90.0 | 40 | 25 | 40 |
| 1233 | 2010-12-05 | 2128 | 42.44 | 71492925.0 | 1.57 | 38.837833 | -122.829833 | 1.59 | 160.0 | 45 | -90.0 | 35 | 20 | 40 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1529 | 2011-03-28 | 732 | 18.81 | 71543515.0 | 1.26 | 38.839000 | -122.825833 | 1.83 | 35.0 | 45 | -100.0 | 25 | 15 | 25 |
| 1530 | 2011-03-29 | 251 | 11.14 | 71543905.0 | 1.38 | 38.838167 | -122.829167 | 1.83 | 90.0 | 55 | -30.0 | 10 | 28 | 30 |
| 1535 | 2011-03-31 | 949 | 8.72 | 71545125.0 | 1.01 | 38.840833 | -122.830500 | 1.75 | 70.0 | 30 | -70.0 | 50 | 23 | 30 |
| 1536 | 2011-03-31 | 949 | 8.72 | 71545125.0 | 1.01 | 38.840833 | -122.830500 | 1.75 | 195.0 | 70 | -140.0 | 25 | 10 | 20 |
| 1537 | 2011-03-31 | 1720 | 8.99 | 71545285.0 | 1.37 | 38.844333 | -122.826500 | 1.84 | 335.0 | 55 | -120.0 | 33 | 20 | 25 |
116 rows × 14 columns
Plot data set
[4]:
def add_scale_bar(ax, x_start, y_start, distance, source_crs, azimuth=90.0, **kwargs):
"""
Parameters
-----------
ax: GeoAxes instance
The axis on which we want to add a scale bar.
x_start: float
The x coordinate of the left end of the scale bar,
given in the axis coordinate system, i.e. from 0 to 1.
y_start: float
The y coordinate of the left end of the scale bar,
given in the axis coordinate system, i.e. from 0 to 1.
distance: float
The distance covered by the scale bar, in km.
source_crs: cartopy.crs
The coordinate system in which the data are written.
"""
from cartopy.geodesic import Geodesic
G = Geodesic()
# default values
kwargs["lw"] = kwargs.get("lw", 2)
kwargs["color"] = kwargs.get("color", "k")
data_coords = ctp.crs.PlateCarree()
# transform the axis coordinates into display coordinates
display = ax.transAxes.transform([x_start, y_start])
# take display coordinates into data coordinates
data = ax.transData.inverted().transform(display)
# take data coordinates into lon/lat
lon_start, lat_start = data_coords.transform_point(data[0], data[1], source_crs)
# get the coordinates of the end of the scale bar
lon_end, lat_end, _ = np.asarray(
G.direct([lon_start, lat_start], azimuth, 1000.0 * distance)
)[0]
if azimuth % 180.0 == 90.0:
lat_end = lat_start
elif azimuth % 180.0 == 0.0:
lon_end == lon_start
ax.plot([lon_start, lon_end], [lat_start, lat_end], transform=data_coords, **kwargs)
ax.text(
(lon_start + lon_end) / 2.0,
(lat_start + lat_end) / 2.0 - 0.001,
"{:.0f}km".format(distance),
transform=data_coords,
ha="center",
va="top",
)
return
def plot_Kaverina(strikes, dips, rakes, ax=None):
n_earthquakes = len(strikes)
# determine the t and p axes
P_axis = np.zeros((n_earthquakes, 2), dtype=np.float32)
T_axis = np.zeros((n_earthquakes, 2), dtype=np.float32)
B_axis = np.zeros((n_earthquakes, 2), dtype=np.float32)
faulting_type = np.zeros(n_earthquakes, dtype=np.int32)
fm_type = []
r2d = 180.0 / np.pi
for t in range(n_earthquakes):
# first, get normal and slip vectors from
# strike, dip, rake
normal, slip = ILSI.utils_stress.normal_slip_vectors(
strikes[t], dips[t], rakes[t]
)
# second, get the t and p vectors
p_axis, t_axis, b_axis = ILSI.utils_stress.p_t_b_axes(normal, slip)
p_bearing, p_plunge = ILSI.utils_stress.get_bearing_plunge(p_axis)
t_bearing, t_plunge = ILSI.utils_stress.get_bearing_plunge(t_axis)
b_bearing, b_plunge = ILSI.utils_stress.get_bearing_plunge(b_axis)
P_axis[t, :] = p_bearing, p_plunge
T_axis[t, :] = t_bearing, t_plunge
B_axis[t, :] = b_bearing, b_plunge
fm_type.append(mecclass(t_plunge, b_plunge, p_plunge))
# get the x, y coordinates for FMC's plots
x_kave, y_kave = kave(T_axis[:, 1], B_axis[:, 1], P_axis[:, 1])
colors_fm = {}
colors_fm["SS"] = "C0"
colors_fm["SS-N"] = "C0"
colors_fm["SS-R"] = "C0"
colors_fm["R"] = "C1"
colors_fm["R-SS"] = "C1"
colors_fm["N"] = "C2"
colors_fm["N-SS"] = "C2"
colors = np.asarray([colors_fm[fm] for fm in fm_type])
if ax is None:
fig = plt.figure("Kaverina_diagram", figsize=(12, 12))
ax = fig.add_subplot(111)
fig = plotFMC.baseplot(10, "", ax=ax)
ax.scatter(x_kave, y_kave, color=colors)
return fig
[5]:
# define region extent
lat_min_box, lat_max_box = 38.83, 38.85
lon_min_box, lon_max_box = -122.84, -122.82
# define inset extent
lat_min_inset, lat_max_inset = 35.00, 42.00
lon_min_inset, lon_max_inset = -125.00, -115.00
cNorm_box = Normalize(
vmin=np.floor(data["depth"].min()), vmax=np.ceil(data["depth"].max())
)
scalar_map_box = ScalarMappable(norm=cNorm_box, cmap=cc.cm.fire_r)
# define projections
data_coords = ctp.crs.PlateCarree()
projection = ctp.crs.Mercator(
central_longitude=(lon_min_box + lon_max_box) / 2.0,
min_latitude=lat_min_box,
max_latitude=lat_max_box,
)
projection_inset = ctp.crs.Mercator(
central_longitude=(lon_min_inset + lon_max_inset) / 2.0,
min_latitude=lat_min_inset,
max_latitude=lat_max_inset,
)
[6]:
# plot the focal mechanisms
fig_fm = plt.figure("earthquakes_Geysers", figsize=(15, 12))
gs = fig_fm.add_gridspec(2, 5)
ax_fm = fig_fm.add_subplot(gs[:, :3], projection=projection)
ax_fm.set_rasterization_zorder(1.5)
ax_fm.set_extent([lon_min_box, lon_max_box, lat_min_box, lat_max_box], crs=data_coords)
ax_fm.set_title("Geysers", pad=15)
# plot inset map
axins = inset_axes(
ax_fm,
width="40%",
height="30%",
loc="upper left",
axes_class=ctp.mpl.geoaxes.GeoAxes,
axes_kwargs=dict(map_projection=projection_inset),
)
axins.set_rasterization_zorder(1.5)
axins.set_extent(
[lon_min_inset, lon_max_inset, lat_min_inset, lat_max_inset], crs=data_coords
)
study_region = geometry.box(
minx=lon_min_box, maxx=lon_max_box, miny=lat_min_box, maxy=lat_max_box
)
axins.add_geometries([study_region], crs=data_coords, edgecolor="k", facecolor="C0")
axins.plot(
(lon_min_box + lon_max_box) / 2.0,
(lat_min_box + lat_max_box) / 2.0,
marker="s",
color="C0",
markersize=10,
markeredgecolor="k",
transform=data_coords,
)
axins.text(
(lon_min_box + lon_max_box) / 2.0 - 0.35,
(lat_min_box + lat_max_box) / 2.0,
"Study Region",
transform=data_coords,
)
# plot San Francisco
axins.plot(
-122.446355,
37.774386,
marker="o",
color="C3",
markersize=10,
markeredgecolor="k",
transform=data_coords,
)
axins.text(-122.15, 37.774386, "San Francisco", transform=data_coords)
for ax in [axins, ax_fm]:
ax.add_feature(
ctp.feature.NaturalEarthFeature(
"cultural",
"admin_1_states_provinces_lines",
"110m",
edgecolor="gray",
facecolor="none",
)
)
ax.add_feature(ctp.feature.BORDERS)
ax.add_feature(ctp.feature.GSHHSFeature(scale="full", levels=[1], zorder=0.49))
axins.add_feature(ctp.feature.LAND)
axins.add_feature(ctp.feature.OCEAN)
add_scale_bar(axins, 0.50, 0.10, 250.0, projection_inset)
# add meridians and latitudes
gl = ax_fm.gridlines(
draw_labels=True, linewidth=1, alpha=0.5, color="k", linestyle="--"
)
gl.right_labels = False
gl.left_labels = True
gl.top_labels = False
gl.bottom_labels = True
for i in range(len(data)):
# fetch focal mechanism
strike, dip, rake = (
data.iloc[i]["strike"],
data.iloc[i]["dip"],
data.iloc[i]["rake"],
)
# determine coordinates in the axis frame of reference
x, y = ax_fm.projection.transform_point(
data.iloc[i]["lon"], data.iloc[i]["lat"], data_coords
)
fc = scalar_map_box.to_rgba(data.iloc[i]["depth"])
# ---------------------------
# uncomment if you are not planning on saving svg figures:
# bb = obb.beach(
# [strike, dip, rake], xy=(x, y), width=20,
# facecolor=fc, linewidth=0.4, axes=ax_fm)
# ---------------------------
# uncomment if you are planning on making svg figures
bb = obb.beach(
[strike, dip, rake],
xy=(x, y),
width=100,
facecolor=fc,
linewidth=0.4,
zorder=1.4,
)
ax_fm.add_collection(bb)
# add colorbar
divider = make_axes_locatable(ax_fm)
cax = divider.append_axes("right", size="5%", pad=0.08, axes_class=plt.Axes)
plt.colorbar(scalar_map_box, cax=cax, label="Depth (km)", orientation="vertical")
# add scale bar
add_scale_bar(ax_fm, 0.05, 0.1, 1.0, projection)
# ---------------------------------------------------
# plot data points in P/T space
ax_pt = fig_fm.add_subplot(gs[0, 3:], projection="stereonet")
# fetch focal mechanisms
strikes, dips, rakes = data["strike"], data["dip"], data["rake"]
n, s = ILSI.utils_stress.normal_slip_vectors(strikes, dips, rakes)
# compute the bearing and plunge of each P/T vector
p_or = np.zeros((len(strikes), 2), dtype=np.float32)
t_or = np.zeros((len(strikes), 2), dtype=np.float32)
for i in range(len(strikes)):
p, t, b = ILSI.utils_stress.p_t_b_axes(n[:, i], s[:, i])
p_or[i, :] = ILSI.utils_stress.get_bearing_plunge(p)
t_or[i, :] = ILSI.utils_stress.get_bearing_plunge(t)
ax_pt.line(t_or[:, 1], t_or[:, 0], marker="v", color="C0")
ax_pt.line(p_or[:, 1], p_or[:, 0], marker="o", color="C3")
ax_pt.line(t_or[0, 1], t_or[0, 0], marker="v", color="C0", label="T-axis")
ax_pt.line(p_or[0, 1], p_or[0, 0], marker="o", color="C3", label="P-axis")
ax_pt.grid(True)
ax_pt.legend(loc="lower left", bbox_to_anchor=(0.80, 0.92))
# ---------------------------------------------------
# plot data points in Kaverina diagram
ax_Kav = fig_fm.add_subplot(gs[1, 3:])
fig_fm = plot_Kaverina(strikes.values, dips.values, rakes.values, ax=ax_Kav)
for i, ax in enumerate([ax_fm, ax_pt, ax_Kav]):
ax.text(
-0.1,
1.05,
f"({string.ascii_lowercase[i]})",
va="top",
fontsize=20,
ha="left",
transform=ax.transAxes,
)
plt.subplots_adjust(wspace=0.55, left=0.1, right=0.98)
ax_pt._polar.set_position(ax_pt.get_position())
/home/eric/miniconda3/envs/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:403: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
warnings.warn("The `map_projection` keyword argument is "
[7]:
if save_figures:
fig_fm.savefig(f'{fig_fm._label}.svg', format='svg', bbox_inches='tight')
Stress tensor inversion
[8]:
# --------------------------------
# prepare data
# --------------------------------
strikes, dips, rakes = data["strike"], data["dip"], data["rake"]
strikes_1, dips_1, rakes_1 = strikes.values, dips.values, rakes.values
rakes_1 = np.float32(rakes_1) % 360.0
n_earthquakes = len(strikes_1)
strikes_2, dips_2, rakes_2 = np.asarray(
list(map(ILSI.utils_stress.aux_plane, strikes_1, dips_1, rakes_1))
).T
[9]:
# --------------------------------
# inversion parameters
# --------------------------------
# if you want to search for the "best" friction coefficient, set friction_coefficient to None
friction_min = 0.1
friction_max = 1.2
friction_step = 0.05
# if you want to fix the friction coefficient to some value, set friction_coefficient equal to that value
# (the paper uses a fixed friction_coefficient = 0.6)
friction_coefficient = 0.6
n_random_selections = 30
n_stress_iter = 10
n_resamplings = 2000
n_averaging = 3
# parallelization option: use n_threads="all" if you want to use all your CPUs to speed up the computation
# use n_threads=X with X any integer number to use a specific number of CPUs
n_threads = "all"
ILSI_kwargs = {}
ILSI_kwargs["max_n_iterations"] = 300
ILSI_kwargs["shear_update_atol"] = 1.0e-5
Tarantola_kwargs = {}
[10]:
# --------------------------------
# initialize output dictionary
# --------------------------------
inversion_output = {}
methods = ["linear", "iterative_constant_shear", "iterative_variable_shear"]
for method in methods:
inversion_output[method] = {}
[11]:
# --------------------------------
# invert the whole dataset
# --------------------------------
print(f"Linear inversion...")
# simple, linear inversion
inversion_output["linear"] = ILSI.ilsi.inversion_one_set(
strikes_1,
dips_1,
rakes_1,
n_random_selections=n_random_selections,
**ILSI_kwargs,
Tarantola_kwargs=Tarantola_kwargs,
variable_shear=False,
return_stats=True,
)
print(f"Iterative constant shear inversion...")
inversion_output["iterative_constant_shear"] = ILSI.ilsi.inversion_one_set_instability(
strikes_1,
dips_1,
rakes_1,
n_random_selections=n_random_selections,
**ILSI_kwargs,
Tarantola_kwargs=Tarantola_kwargs,
friction_min=friction_min,
friction_max=friction_max,
friction_step=friction_step,
friction_coefficient=friction_coefficient,
n_stress_iter=n_stress_iter,
n_averaging=n_averaging,
variable_shear=False,
signed_instability=False,
return_stats=True,
)
if friction_coefficient is None:
print(
f'Inverted friction: {inversion_output["iterative_constant_shear"]["friction_coefficient"]:.2f}'
)
print(f"Iterative variable shear inversion...")
inversion_output["iterative_variable_shear"] = ILSI.ilsi.inversion_one_set_instability(
strikes_1,
dips_1,
rakes_1,
n_random_selections=n_random_selections,
**ILSI_kwargs,
Tarantola_kwargs=Tarantola_kwargs,
friction_min=friction_min,
friction_max=friction_max,
friction_step=friction_step,
friction_coefficient=friction_coefficient,
n_stress_iter=n_stress_iter,
n_averaging=n_averaging,
variable_shear=True,
signed_instability=False,
return_stats=True,
plot=False,
)
if friction_coefficient is None:
print(
f'Inverted friction: {inversion_output["iterative_variable_shear"]["friction"]:.2f}'
)
for method in methods:
R = ILSI.utils_stress.R_(inversion_output[method]["principal_stresses"])
I, fp_strikes, fp_dips, fp_rakes = ILSI.ilsi.compute_instability_parameter(
inversion_output[method]["principal_directions"],
R,
inversion_output[method]["friction_coefficient"] if "friction_coefficient" in inversion_output[method] else 0.6,
strikes_1,
dips_1,
rakes_1,
strikes_2,
dips_2,
rakes_2,
return_fault_planes=True,
)
inversion_output[method]["misfit"] = np.mean(
ILSI.utils_stress.mean_angular_residual(
inversion_output[method]["stress_tensor"], fp_strikes, fp_dips, fp_rakes
)
)
Linear inversion...
Iterative constant shear inversion...
-------- 1/3 ----------
Initial shape ratio: 0.43
-------- 2/3 ----------
Initial shape ratio: 0.43
-------- 3/3 ----------
Initial shape ratio: 0.43
Final results:
Stress tensor:
[[ 0.1799634 0.26173306 -0.19272774]
[ 0.26173306 0.47864988 0.21612422]
[-0.19272774 0.21612422 -0.65861326]]
Shape ratio: 0.63
Iterative variable shear inversion...
-------- 1/3 ----------
Initial shape ratio: 0.47
-------- 2/3 ----------
Initial shape ratio: 0.47
-------- 3/3 ----------
Initial shape ratio: 0.48
Final results:
Stress tensor:
[[ 0.28717783 0.22700854 -0.24951024]
[ 0.22700854 0.5183292 0.30055955]
[-0.24951024 0.30055955 -0.805507 ]]
Shape ratio: 0.75
[12]:
# --------------------------------
# bootstrap the dataset to infer uncertainties
# --------------------------------
print(f"Linear inversion (bootstrapping)...")
t1 = give_time()
# simple, linear inversion
bootstrap = ILSI.ilsi.inversion_bootstrap(
strikes_1,
dips_1,
rakes_1,
n_resamplings=n_resamplings,
Tarantola_kwargs=Tarantola_kwargs,
variable_shear=False,
)
for field in bootstrap:
inversion_output["linear"][field] = bootstrap[field]
t2 = give_time()
print(f"Done in {t2-t1:.2f}sec!")
t1 = give_time()
print(f"Iterative constant shear inversion (bootstrapping)...")
bootstrap = ILSI.ilsi.inversion_bootstrap_instability(
inversion_output["iterative_constant_shear"]["principal_directions"],
ILSI.utils_stress.R_(
inversion_output["iterative_constant_shear"]["principal_stresses"]
),
strikes_1,
dips_1,
rakes_1,
inversion_output["iterative_constant_shear"]["friction_coefficient"],
**ILSI_kwargs,
Tarantola_kwargs=Tarantola_kwargs,
n_stress_iter=n_stress_iter,
n_resamplings=n_resamplings,
n_threads=n_threads,
variable_shear=False,
signed_instability=False,
)
for field in bootstrap:
inversion_output["iterative_constant_shear"][field] = bootstrap[field]
t2 = give_time()
print(f"Done in {t2-t1:.2f}sec!")
t1 = give_time()
print(f"Iterative variable shear inversion (bootstrapping)...")
bootstrap = ILSI.ilsi.inversion_bootstrap_instability(
inversion_output["iterative_variable_shear"]["principal_directions"],
ILSI.utils_stress.R_(
inversion_output["iterative_variable_shear"]["principal_stresses"]
),
strikes_1,
dips_1,
rakes_1,
inversion_output["iterative_variable_shear"]["friction_coefficient"],
**ILSI_kwargs,
Tarantola_kwargs=Tarantola_kwargs,
n_stress_iter=n_stress_iter,
n_resamplings=n_resamplings,
n_threads=n_threads,
variable_shear=True,
signed_instability=False,
)
for field in bootstrap:
inversion_output["iterative_variable_shear"][field] = bootstrap[field]
t2 = give_time()
print(f"Done in {t2-t1:.2f}sec!")
Linear inversion (bootstrapping)...
---------- Bootstrapping 1/2000 ----------
---------- Bootstrapping 101/2000 ----------
---------- Bootstrapping 201/2000 ----------
---------- Bootstrapping 301/2000 ----------
---------- Bootstrapping 401/2000 ----------
---------- Bootstrapping 501/2000 ----------
---------- Bootstrapping 601/2000 ----------
---------- Bootstrapping 701/2000 ----------
---------- Bootstrapping 801/2000 ----------
---------- Bootstrapping 901/2000 ----------
---------- Bootstrapping 1001/2000 ----------
---------- Bootstrapping 1101/2000 ----------
---------- Bootstrapping 1201/2000 ----------
---------- Bootstrapping 1301/2000 ----------
---------- Bootstrapping 1401/2000 ----------
---------- Bootstrapping 1501/2000 ----------
---------- Bootstrapping 1601/2000 ----------
---------- Bootstrapping 1701/2000 ----------
---------- Bootstrapping 1801/2000 ----------
---------- Bootstrapping 1901/2000 ----------
Done in 9.43sec!
Iterative constant shear inversion (bootstrapping)...
Done in 9.08sec!
Iterative variable shear inversion (bootstrapping)...
Done in 9.83sec!
[13]:
inversion_output["strikes"] = np.stack((strikes_1, strikes_2), axis=1)
inversion_output["dips"] = np.stack((dips_1, dips_2), axis=1)
inversion_output["rakes"] = np.stack((rakes_1, rakes_2), axis=1)
[14]:
def plot_inverted_stress_tensors(
inversion,
axes=None,
figtitle="",
figname="inverted_stress_tensor",
methods=["linear", "iterative_constant_shear", "iterative_variable_shear"],
figsize=(15, 11),
plot_PT=False,
error="bootstrap",
**kwargs,
):
"""
Plots the inverted stress tensors from a given stress inversion analysis.
Parameters
----------
inversion : dict
Dictionary containing the results of the stress tensor inversion.
It should include principal stresses, stress tensors, and bootstrapped
principal directions for different inversion methods.
axes : list, optional
List of axes to plot on. If None, new axes will be created.
figtitle : str, optional
Title of the figure. Default is an empty string.
figname : str, optional
Name of the figure window. Default is "inverted_stress_tensor".
methods : list of str, optional
List of stress inversion methods to plot.
Default is ["linear", "iterative_constant_shear", "iterative_variable_shear"].
figsize : tuple of (float, float), optional
Figure size in inches. Default is (15, 11).
plot_PT : bool, optional
If True, plots P- and T-axes from fault slip data. Default is False.
error : {'bootstrap', 'covariance'}, optional
Specifies the error estimation method to use:
- 'bootstrap': Uses bootstrapped principal directions.
- 'covariance': Uses covariance matrix to sample stress tensors.
Default is 'bootstrap'.
**kwargs : dict, optional
Additional parameters for histogram plotting, including:
- smoothing_sig : float, optional
Smoothing parameter for density estimation. Default is 1.
- nbins : int, optional
Number of bins for histograms. Default is 200.
- return_count : bool, optional
If True, returns count data for histograms. Default is True.
- confidence_intervals : list of float, optional
List of confidence intervals (in percentages). Default is [95.0].
Returns
-------
fig : matplotlib.figure.Figure
The generated matplotlib figure containing the stress tensor plots.
"""
hist_kwargs = {}
hist_kwargs["smoothing_sig"] = kwargs.get("smoothing_sig", 1)
hist_kwargs["nbins"] = kwargs.get("nbins", 200)
hist_kwargs["return_count"] = kwargs.get("return_count", True)
hist_kwargs["confidence_intervals"] = kwargs.get("confidence_intervals", [95.0])
markers = ["o", "s", "v"]
n_resamplings = inversion[methods[0]]["boot_principal_directions"].shape[0]
strikes = inversion["strikes"]
dips = inversion["dips"]
rakes = inversion["rakes"]
fig = plt.figure(figname, figsize=figsize)
fig.suptitle(figtitle)
gs = fig.add_gridspec(
nrows=3,
ncols=3,
top=0.88,
bottom=0.08,
left=0.15,
right=0.95,
hspace=0.4,
wspace=0.7,
)
axes = []
ax1 = fig.add_subplot(gs[0, 0], projection="stereonet")
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
for j, method in enumerate(methods):
if error == "bootstrap":
replica_principal_directions = inversion[method]["boot_principal_directions"]
replica_principal_stresses = inversion[method]["boot_principal_stresses"]
elif error == "covariance":
st = inversion[method]["stress_tensor"]
sigma_mean = st[np.triu_indices(3)][:-1]
sigma_cov = inversion[method]["C_m_posterior"]
# correct symmetry for numerical errors
sigma_cov = (sigma_cov + sigma_cov.T) / 2.
random_stress_tensors = np.random.multivariate_normal(
sigma_mean, sigma_cov, size=n_resamplings
)
random_stress_tensors = np.stack(
[
random_stress_tensors[:, 0],
random_stress_tensors[:, 1],
random_stress_tensors[:, 2],
random_stress_tensors[:, 1],
random_stress_tensors[:, 3],
random_stress_tensors[:, 4],
random_stress_tensors[:, 2],
random_stress_tensors[:, 4],
-random_stress_tensors[:, 0] - random_stress_tensors[:, 3],
],
axis=1
).reshape(n_resamplings, 3, 3)
replica_principal_directions = np.zeros_like(random_stress_tensors)
replica_principal_stresses = np.zeros((n_resamplings, 3), dtype=np.float32)
for n in range(n_resamplings):
replica_principal_stresses[n, :], replica_principal_directions[n, ...] = ILSI.utils_stress.stress_tensor_eigendecomposition(
random_stress_tensors[n, ...]
)
R = ILSI.utils_stress.R_(inversion[method]["principal_stresses"])
A_phi = ILSI.utils_stress.A_phi_(
inversion[method]["principal_stresses"],
inversion[method]["principal_directions"],
)
for k in range(3):
if k == 0:
label = "{}:\n".format(method.replace("_", " ").capitalize())
for k2 in range(3):
az, pl = ILSI.utils_stress.get_bearing_plunge(
inversion[method]["principal_directions"][:, k2]
)
label += (
r"$\sigma_{{{:d}}}$: {:.1f}"
"\u00b0"
"|{:.1f}"
"\u00b0"
", ".format(k2 + 1, az, pl)
)
label += (
f" R={R:.2f}"
r", $A_{\phi}=$"
f"{A_phi:.2f}, $\\beta $={inversion[method]['misfit']:.1f}"
"\u00b0"
)
else:
label = ""
az, pl = ILSI.utils_stress.get_bearing_plunge(
inversion[method]["principal_directions"][:, k]
)
ax1.line(
pl,
az,
marker=markers[k],
markeredgecolor=_colors_[j],
markeredgewidth=2,
markerfacecolor="none",
markersize=[15, 15, 15, 15][j],
label=label,
zorder=2,
)
boot_pd_stereo = np.zeros((n_resamplings, 2), dtype=np.float32)
for b in range(n_resamplings):
boot_pd_stereo[b, :] = ILSI.utils_stress.get_bearing_plunge(
replica_principal_directions[b, :, k]
)
count, lons_g, lats_g, levels = ILSI.utils_stress.get_CI_levels(
boot_pd_stereo[:, 0], boot_pd_stereo[:, 1], **hist_kwargs
)
ax1.contour(
lons_g,
lats_g,
count,
levels=levels,
vmin=0.0,
linestyles=["solid", "dashed", "dashdot"][k],
linewidths=0.75,
colors=_colors_[j],
zorder=1.1,
)
if plot_PT:
# get fault noral and slip vector from strikes/dips/rakes
n, s = ILSI.utils_stress.normal_slip_vectors(strikes[:, 0], dips[:, 0], rakes[:, 0])
# compute the bearing and plunge of each P/T vector
p_or = np.zeros((strikes.shape[0], 2), dtype=np.float32)
t_or = np.zeros((strikes.shape[0], 2), dtype=np.float32)
for i in range(strikes.shape[0]):
p, t, b = ILSI.utils_stress.p_t_b_axes(n[:, i], s[:, i])
p_or[i, :] = ILSI.utils_stress.get_bearing_plunge(p)
t_or[i, :] = ILSI.utils_stress.get_bearing_plunge(t)
ax1.line(t_or[:, 1], t_or[:, 0], zorder=0.5, marker="v", color="grey")
ax1.line(p_or[:, 1], p_or[:, 0], zorder=0.5, marker="o", color="k")
if j == len(methods) - 1:
ax1.line(t_or[0, 1], t_or[0, 0], zorder=0.5, marker="v", color="grey", label="T-axis")
ax1.line(p_or[0, 1], p_or[0, 0], zorder=0.5, marker="o", color="k", label="P-axis")
axes.append(ax1)
Rs = np.zeros(n_resamplings, dtype=np.float32)
for b in range(n_resamplings):
Rs[b] = ILSI.utils_stress.R_(
replica_principal_stresses[b, :]
)
ax2.hist(
Rs, range=(0.0, 1.0), bins=20, lw=2.5, color=_colors_[j], histtype="step"
)
ax2.set_xlabel("Shape Ratio")
ax2.set_ylabel("Count")
# plot stress magnitudes
for i, method in enumerate(methods):
pd = inversion[method]["principal_directions"]
R = ILSI.utils_stress.R_(inversion[method]["principal_stresses"])
I, fp_strikes, fp_dips, fp_rakes = ILSI.ilsi.compute_instability_parameter(
pd,
R,
inversion_output[method]["friction_coefficient"] if "friction_coefficient" in inversion_output[method] else 0.6,
strikes[:, 0],
dips[:, 0],
rakes[:, 0],
strikes[:, 1],
dips[:, 1],
rakes[:, 1],
return_fault_planes=True,
)
normals, slips = ILSI.utils_stress.normal_slip_vectors(
fp_strikes, fp_dips, fp_rakes
)
_, _, shear = ILSI.utils_stress.compute_traction(
inversion[method]["stress_tensor"], normals.T
)
Ts_mag = np.sqrt(np.sum(shear**2, axis=-1))
sig_Ts = np.std(Ts_mag)
ax3.hist(
Ts_mag,
range=(0.0, 1.0),
bins=20,
color=_colors_[i],
alpha=0.5,
label=method.replace("_", " ").capitalize(),
)
ax3.set_ylabel("Count")
ax3.set_xlabel("Shear stress $\\vert \\tau \\vert$")
ax3.legend(bbox_to_anchor=(-0.25, 1.05), loc="lower left")
ax1.legend(loc="upper left", bbox_to_anchor=(-0.5, -0.30))
for i, ax in enumerate([ax1, ax2, ax3]):
ax.text(
-0.25,
1.2,
f"({string.ascii_lowercase[i]})",
va="top",
fontsize=20,
ha="left",
transform=ax.transAxes,
)
plt.subplots_adjust(
top=0.88, bottom=0.11, left=0.05, right=0.95, hspace=0.4, wspace=0.4
)
ax1._polar.set_position(ax1.get_position())
return fig
[27]:
# plot results
kwargs = {}
# -----------------------------------------------------------------
# check out the confidence intervals given by bootstrapping
# -----------------------------------------------------------------
kwargs["smoothing_sig"] = 5 # control the smoothness of the confidence intervals
kwargs["confidence_intervals"] = [95.0]
fig_SI = plot_inverted_stress_tensors(
inversion_output, error="bootstrap", figtitle="Geysers data set", **kwargs
)
# -----------------------------------------------------------------
# check out the confidence intervals given by the Tarantola-Valette solution
# -----------------------------------------------------------------
# kwargs["smoothing_sig"] = 5 # control the smoothness of the confidence intervals
# kwargs["confidence_intervals"] = [66.0]
# fig_SI = plot_inverted_stress_tensors(
# inversion_output, error="covariance", figtitle="Geysers data set", **kwargs
# )
[16]:
# save the figures
if save_figures:
for fig in [fig_SI]:
fig.savefig(f'{fig._label}.png', format='png', bbox_inches='tight')
fig.savefig(f'{fig._label}.svg', format='svg', bbox_inches='tight')
[17]:
def print_CI(inversion_output, CI_level=95.0, **kwargs):
hist_kwargs = {}
hist_kwargs["smoothing_sig"] = kwargs.get("smoothing_sig", 1)
hist_kwargs["nbins"] = kwargs.get("nbins", 200)
hist_kwargs["return_count"] = kwargs.get("return_count", True)
hist_kwargs["confidence_intervals"] = [CI_level]
methods = ["linear", "iterative_constant_shear", "iterative_variable_shear"]
n_resamplings = inversion_output["linear"]["boot_principal_directions"].shape[0]
for j in [0, 1, 2]:
method = methods[j]
print(f"---------- {method} ---------")
R = ILSI.utils_stress.R_(inversion_output[method]["principal_stresses"])
Rs = np.zeros(n_resamplings, dtype=np.float32)
for b in range(n_resamplings):
Rs[b] = ILSI.utils_stress.R_(
inversion_output[method]["boot_principal_stresses"][b, :]
)
R_minus = np.percentile(Rs, (100.0 - CI_level) / 2.0)
R_plus = np.percentile(Rs, CI_level + (100.0 - CI_level) / 2.0)
print(
f"\tShape ratio R = {R:.2f}, CI {CI_level:.0f} = ({R_minus:.2f}, {R_plus:.2f})"
)
for k in range(3):
az, pl = ILSI.utils_stress.get_bearing_plunge(
inversion_output[method]["principal_directions"][:, k]
)
boot_pd_stereo = np.zeros((n_resamplings, 2), dtype=np.float32)
for b in range(n_resamplings):
boot_pd_stereo[b, :] = ILSI.utils_stress.get_bearing_plunge(
inversion_output[method]["boot_principal_directions"][b, :, k]
)
count, lons_g, lats_g, levels = ILSI.utils_stress.get_CI_levels(
boot_pd_stereo[:, 0], boot_pd_stereo[:, 1], **hist_kwargs
)
plunges, bearings = mpl.geographic2plunge_bearing(lons_g, lats_g)
CR = np.where(count.flatten() >= levels[0])[0]
CR_plunges = plunges.flatten()[CR]
CR_bearings = bearings.flatten()[CR]
# correct for possible "spread" of the CI at the other
# end of the line (+/- 180deg) for low plunges
CR_bearings[CR_plunges < 15.0] = CR_bearings[CR_plunges < 15.0] % 180.0
if az > 180.0:
CR_bearings[CR_plunges < 15.0] += 180.0
print(f"\tPrincipal stress sigma {k+1}:")
print(
f"\t\tPlunge = {pl:.2f}"
"\u00b0"
f", CI {CI_level:.0f} = ({CR_plunges.min():.2f}"
"\u00b0"
f", {CR_plunges.max():.2f}"
"\u00b0"
f")"
)
print(
f"\t\tBearing = {az:.2f}"
"\u00b0"
f", CI {CI_level:.0f} = ({CR_bearings.min():.2f}"
"\u00b0"
f", {CR_bearings.max():.2f}"
"\u00b0"
f")"
)
print(
f"\t\t({az:.1f}"
"\u00b0"
f", {pl:.1f}"
"\u00b0"
f") ({CR_bearings.min():.1f}"
"\u00b0"
f"-{CR_bearings.max():.1f}"
"\u00b0"
f", {CR_plunges.min():.1f}"
"\u00b0"
f"-{CR_plunges.max():.1f}"
"\u00b0"
f")"
)
[18]:
kwargs = {}
kwargs["smoothing_sig"] = 5 # control the smoothness of the confidence intervals
kwargs["confidence_intervals"] = [95.0]
print_CI(inversion_output, **kwargs)
---------- linear ---------
Shape ratio R = 0.43, CI 95 = (0.30, 0.57)
Principal stress sigma 1:
Plunge = 66.99°, CI 95 = (51.14°, 83.61°)
Bearing = 216.63°, CI 95 = (181.50°, 257.14°)
(216.6°, 67.0°) (181.5°-257.1°, 51.1°-83.6°)
Principal stress sigma 2:
Plunge = 22.39°, CI 95 = (8.73°, 38.84°)
Bearing = 22.54°, CI 95 = (10.77°, 38.71°)
(22.5°, 22.4°) (10.8°-38.7°, 8.7°-38.8°)
Principal stress sigma 3:
Plunge = 5.05°, CI 95 = (0.34°, 16.53°)
Bearing = 114.63°, CI 95 = (101.25°, 130.06°)
(114.6°, 5.1°) (101.3°-130.1°, 0.3°-16.5°)
---------- iterative_constant_shear ---------
Shape ratio R = 0.63, CI 95 = (0.44, 0.76)
Principal stress sigma 1:
Plunge = 70.46°, CI 95 = (53.75°, 85.55°)
Bearing = 220.66°, CI 95 = (177.54°, 267.26°)
(220.7°, 70.5°) (177.5°-267.3°, 53.8°-85.5°)
Principal stress sigma 2:
Plunge = 19.00°, CI 95 = (5.99°, 35.80°)
Bearing = 26.64°, CI 95 = (9.74°, 44.96°)
(26.6°, 19.0°) (9.7°-45.0°, 6.0°-35.8°)
Principal stress sigma 3:
Plunge = 4.39°, CI 95 = (0.32°, 17.56°)
Bearing = 118.16°, CI 95 = (101.26°, 134.59°)
(118.2°, 4.4°) (101.3°-134.6°, 0.3°-17.6°)
---------- iterative_variable_shear ---------
Shape ratio R = 0.75, CI 95 = (0.49, 0.82)
Principal stress sigma 1:
Plunge = 70.78°, CI 95 = (58.96°, 89.36°)
Bearing = 224.56°, CI 95 = (135.00°, 315.04°)
(224.6°, 70.8°) (135.0°-315.0°, 59.0°-89.4°)
Principal stress sigma 2:
Plunge = 18.22°, CI 95 = (0.09°, 31.60°)
Bearing = 25.36°, CI 95 = (3.31°, 42.89°)
(25.4°, 18.2°) (3.3°-42.9°, 0.1°-31.6°)
Principal stress sigma 3:
Plunge = 5.90°, CI 95 = (0.34°, 21.60°)
Bearing = 117.31°, CI 95 = (92.28°, 131.02°)
(117.3°, 5.9°) (92.3°-131.0°, 0.3°-21.6°)
Bonus: Uncertainty analysis beside bootstrapping: The resolution and posterior covariance operators
The previous figure shows bootstrapping uncertainties. Since the Tarantola-Valette inverse method uses the Bayesian framework, among the byproducts of the inversion are the posterior covariance operators \(\hat{C}_D\) and \(\hat{C}_M\). Furthermore, a generic concept to gauge the quality of an inverse method is the resolution operator \(R\). Here, we simply demonstrate the tools implemented in ILSI to manipulate these operators. These concepts are discussed in detail in:
Lundstern, Jens-Erik, Eric Beaucé, and Orlando J. Teran. “The importance of nodal plane orientation diversity for earthquake focal mechanism stress inversions.” Geological Society, London (2024),
and also illustrated in:
In general, bootstrapping will provide more realistic uncertainties but, in some cases where the input data set is small or relatively large but with little focal mechanism diversity, the posterior covariance and resolution operators may provide better uncertainties.
Note: The 5-element \(m\) vector is \((\sigma_{11}, \sigma_{12}, \sigma_{13}, \sigma_{22}, \sigma_{23})\). We work in the (N=north, W=west, U=up) coordinate system, therefore \(m\) is: \((\sigma_{NN}, \sigma_{NW}, \sigma_{NU}, \sigma_{WW}, \sigma_{WU})\)
[19]:
def plot_resolution_covariance_corrcoef(
R,
C_m_post,
C_m_post_normalized,
figname="resolution_covariance",
figtitle=r"Resolution operator $R = LG$, posterior covariance $\hat{C}_M$"
r" and normalized posterior covariance $\rho_M$",
figsize=(13, 10),
cmap="magma"
):
"""Plot resolution operator, posterior covariance, and normalized posterior covariance.
Parameters:
-----------
R : array_like
Resolution operator.
C_m_post : array_like
Posterior covariance.
C_m_post_normalized : array_like
Normalized posterior covariance.
figname : str, optional
Figure name (default is 'resolution_covariance').
figtitle : str, optional
Figure title (default is 'Resolution operator R = LG, posterior
covariance C_M, and normalized posterior covariance ρ_M').
figsize : tuple, optional
Figure size (default is (13, 10)).
cmap : str or Colormap, optional
Colormap for plotting (default is 'viridis').
"""
fig, axes = plt.subplots(num=figname, ncols=2, nrows=2, figsize=figsize)
fig.suptitle(figtitle)
plt.subplots_adjust(top=0.90, bottom=0.08, hspace=0.30)
axes[0, 0].set_title(r"Resolution operator, $R = LG$")
pc0 = axes[0, 0].pcolormesh(R, cmap=cmap, rasterized=True,)
plt.colorbar(pc0, label="Resolution operator")
diff = np.abs(np.identity(5) - R)
axes[0, 1].set_title(r"$\vert I - R \vert$")
pc1 = axes[0, 1].pcolormesh(diff, cmap=cmap, rasterized=True, vmin=0.,)
plt.colorbar(pc1, label="Deviation from identity")
axes[1, 0].set_title(r"Posterior covariance operator, $\hat{C}_M$")
pc2 = axes[1, 0].pcolormesh(C_m_post, cmap=cmap, rasterized=True)
plt.colorbar(pc2, label="Covariance")
vmin = min(C_m_post_normalized[C_m_post_normalized < 0.999].min(), -C_m_post_normalized[C_m_post_normalized < 0.999].max())+0.01
vmax = -vmin
# vmin, vmax = -1., +1.
cm = plt.cm.get_cmap("coolwarm")
cm.set_over("w")
axes[1, 1].set_title(r"Normalized post. covariance op., $\rho_M$")
pc3 = axes[1, 1].pcolormesh(
C_m_post_normalized, rasterized=True, vmin=vmin, vmax=vmax, cmap=cm
)
plt.colorbar(pc3, label="Normalized covariance")
tickpos = [0.5, 1.5, 2.5, 3.5, 4.5]
ticklabels = [
r"$\sigma_{11}$",
r"$\sigma_{12}$",
r"$\sigma_{13}$",
r"$\sigma_{22}$",
r"$\sigma_{23}$"
]
for i, ax in enumerate(axes.flatten()):
ax.set_xlabel("True model parameter")
ax.set_ylabel("Inverted model parameter")
ax.set_xticks(tickpos)
ax.set_xticklabels(ticklabels)
ax.set_yticks(tickpos)
ax.set_yticklabels(ticklabels)
ax.text(
-0.1,
1.05,
f"({string.ascii_lowercase[i]})",
transform=ax.transAxes,
size=20,
)
return fig
[20]:
fig_resolutions = []
for method in ["linear", "iterative_constant_shear", "iterative_variable_shear"]:
C_m_post = inversion_output[method]["C_m_posterior"]
C_m_post_normalized = ( C_m_post / np.sqrt(np.diag(C_m_post))[:, None]) / (np.sqrt(np.diag(C_m_post))[None, :])
fig_resolutions.append(plot_resolution_covariance_corrcoef(
inversion_output[method]["resolution_operator"],
inversion_output[method]["C_m_posterior"],
C_m_post_normalized,
figname=f"resolution_covariance_{method}",
))
fig_resolutions[-1].suptitle(
fig_resolutions[-1].texts[0].get_text() + f"\n{method.replace('_', ' ').capitalize()}", y=1.01
)
Some extra figures
[21]:
def plot_Mohr(inversion, plot_density=False):
import synthetic_dataset as syndata
strikes = inversion["strikes"]
dips = inversion["dips"]
rakes = inversion["rakes"]
fig, axes = plt.subplots(
num=f"Mohr_space_Geysers", figsize=(15, 5), nrows=1, ncols=3
)
fig.suptitle(f"The Geysers geothermal field")
for i, method in enumerate(
["linear", "iterative_constant_shear", "iterative_variable_shear"]
):
axes[i].set_title(method.replace("_", " ").capitalize())
st = inversion[method]["stress_tensor"]
p_sig, p_dir = ILSI.utils_stress.stress_tensor_eigendecomposition(st)
I, s, d, r = ILSI.ilsi.compute_instability_parameter(
p_dir,
ILSI.utils_stress.R_(p_sig),
0.6,
strikes[:, 0],
dips[:, 0],
rakes[:, 0],
strikes[:, 1],
dips[:, 1],
rakes[:, 1],
return_fault_planes=True,
)
fig = syndata.plot_dataset_Mohr(st, s, d, plot_density=plot_density, ax=axes[i])
fig = syndata.plot_dataset_Mohr(st, s, d, plot_density=plot_density, ax=axes[i])
fig = syndata.plot_dataset_Mohr(st, s, d, plot_density=plot_density, ax=axes[i])
for i, ax in enumerate(axes.flatten()):
ax.text(
-0.1,
1.05,
f"({string.ascii_lowercase[i]})",
transform=ax.transAxes,
size=20,
)
fig.tight_layout()
return fig
[24]:
fig_Mohr = plot_Mohr(inversion_output, plot_density=True)
[25]:
def plot_instabilities(inversion, mu=0.60, Imin=0.80):
strikes = inversion["strikes"]
dips = inversion["dips"]
rakes = inversion["rakes"]
fig, ax = plt.subplots(num="instabilities", figsize=(5, 5), nrows=1, ncols=1)
for i, method in enumerate(
["linear", "iterative_constant_shear", "iterative_variable_shear"]
):
stress_tensor = inversion[method]["stress_tensor"]
p_sig, p_dir = ILSI.utils_stress.stress_tensor_eigendecomposition(stress_tensor)
I, s, d, r = ILSI.ilsi.compute_instability_parameter(
p_dir,
ILSI.utils_stress.R_(p_sig),
mu,
strikes[:, 0],
dips[:, 0],
rakes[:, 0],
strikes[:, 1],
dips[:, 1],
rakes[:, 1],
return_fault_planes=True,
)
ax.hist(
np.max(I, axis=-1),
bins=20,
range=(Imin, 1.0),
color=_colors_[i],
label=method.replace("_", " ").capitalize(),
alpha=1,
histtype="step",
linewidth=3.0,
)
ax.set_xlabel(r"Instability $I$")
ax.set_ylabel("Count")
ax.legend(loc="upper left")
return fig
[26]:
fig_inst = plot_instabilities(inversion_output, Imin = 0.)
ax = fig_inst.get_axes()[0]
ax.set_title('Geysers')
fig_inst.savefig(fig_inst._label+'_Geysers.svg', format='svg', bbox_inches='tight')
