Skip to content

DOC: add "strain and stress calculation" in Tutorial, and fix some typos #16

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Apr 18, 2025
Merged

DOC: add "strain and stress calculation" in Tutorial, and fix some typos #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5fa0d36
DOC: add `grt.b2a` example in Tutorial
Dengda98 Apr 18, 2025
351d0d9
DOC: add some formulas for synthetic seismograms in Tutorial
Dengda98 Apr 18, 2025
9627429
DOC: remove some typo
Dengda98 Apr 18, 2025
534ce7d
DOC: add `docs/` reference
Dengda98 Apr 18, 2025
1cc4535
DOC: add "strain and stress calculation" in Tutorial
Dengda98 Apr 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions docs/source/Tutorial/dynamic/gfunc.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,22 @@ Python中计算动态格林函数的主函数为 :func:`compute_grn() <pygrt.pym
:start-after: BEGIN GRN
:end-before: END GRN

不同震源深度、接收点深度和震中距的格林函数会在 :rst:dir:`GRN/milrow_{depsrc}_{deprcv}_{dist}/`` 路径下，使用SAC格式保存。
不同震源深度、接收点深度和震中距的格林函数会在 :rst:dir:`GRN/milrow_{depsrc}_{deprcv}_{dist}/` 路径下，使用SAC格式保存。

一些基本信息（包括源点和场点的物性参数）保存在SAC头段变量中，其中 :c:var:`t0` 和 :c:var:`t1` 分别代表初至P波和初至S波的到时，由于程序中使用0作为参考时间，故其等价于走时。

如果你没有安装SAC软件，可以使用Python的ObsPy库读取数据，或者使用 :command:`grt.b2a` 工具临时将SAC格式文件转为如下的文本文件：

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN grt.b2a
:end-before: END grt.b2a

输出的文本文件如下，两列分别为时间点和幅值。这种输出仅保留波形信息，缺失SAC文件中的头段变量。

.. literalinclude:: run/HFZ_head
:language: text

.. tab:: Python

.. literalinclude:: run/run.py
Expand Down Expand Up @@ -95,7 +107,7 @@ Python中计算动态格林函数的主函数为 :func:`compute_grn() <pygrt.pym
积分形式分类
--------------

通过在面谐矢量坐标系中建立波函数进行公式推导，最终理论地震图的三分量频谱 :math:`W_m(\omega), Q_m(\omega), V_m(\omega)` (分别为垂向，径向，切向)可以表达为:
通过在面谐矢量坐标系中建立波函数进行公式推导，最终格林函数的三分量频谱 :math:`W_m(\omega), Q_m(\omega), V_m(\omega)` (分别为垂向，径向，切向)可以表达为:

.. math::

Expand Down
11 changes: 11 additions & 0 deletions docs/source/Tutorial/dynamic/run/HFZ_head
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
0.0000000e+00 -6.9433704e-06
2.0000000e-02 8.6446698e-06
3.9999999e-02 -8.0442423e-06
5.9999999e-02 9.3830195e-06
7.9999998e-02 -9.2141536e-06
9.9999994e-02 9.3351609e-06
1.2000000e-01 -9.4855750e-06
1.4000000e-01 1.3793133e-05
1.6000000e-01 -8.3573850e-06
1.7999999e-01 1.0849539e-05
...
6 changes: 6 additions & 0 deletions docs/source/Tutorial/dynamic/run/run.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,12 @@
grt -Mmilrow -D2/0 -N500/0.02 -OGRN -R5,8,10
# END GRN

# BEGIN grt.b2a
grt.b2a GRN/milrow_2_0_10/HFZ.sac > HFZ
# END grt.b2a
head -n 10 HFZ > HFZ_head
echo "..." >> HFZ_head

# BEGIN SYN EXP
# 合成结果在 syn_exp/ 目录下，以SAC格式保存。
grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -Osyn_exp
Expand Down
9 changes: 9 additions & 0 deletions docs/source/Tutorial/dynamic/run_upar/milrow
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
0.2 3.4 1.7 2.3 9e10 9e10
0.6 3.7 1.9 2.4 9e10 9e10
0.5 4.2 2.1 2.4 9e10 9e10
0.5 4.6 2.3 2.5 9e10 9e10
0.7 4.9 2.8 2.6 9e10 9e10
0.5 5.1 2.9 2.7 9e10 9e10
6.0 5.9 3.3 2.7 9e10 9e10
1. 6.9 4.0 2.8 9e10 9e10
0.0 8.2 4.7 3.2 9e10 9e10
173 changes: 173 additions & 0 deletions docs/source/Tutorial/dynamic/run_upar/run.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
# BEGIN GRN
import numpy as np
import pygrt

modarr = np.loadtxt("milrow")

pymod = pygrt.PyModel1D(modarr, depsrc=2.0, deprcv=0.0)

# 传入calc_upar=True计算空间导数
stgrn = pymod.compute_grn(distarr=[10], nt=500, dt=0.02, calc_upar=True)[0]
print(stgrn.__str__(extended=True))
# 45 Trace(s) in Stream:
# .SYN..EXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..VFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..HFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..SSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zEXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rEXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zVFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rVFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zHFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rHFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zSSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rSSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..EXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..VFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..HFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..SSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zEXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rEXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zVFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rVFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zHFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rHFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zSSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rSSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..HFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..SST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zHFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rHFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zSST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rSST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# END GRN



# BEGIN SYN DC
# 传入calc_upar=True计算空间导数
stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30, calc_upar=True)
print(stsyn)
# 12 Trace(s) in Stream:
# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..rDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..tDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..tDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..tDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# END SYN DC


# BEGIN ZNE
# 传入ZNE=True可返回ZNE分量
stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30, calc_upar=True, ZNE=True)
print(stsyn)
# 12 Trace(s) in Stream:
# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..DCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..zDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..nDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..nDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..nDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..eDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..eDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..eDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# END ZNE

from obspy import *
import matplotlib.pyplot as plt

def plot6(st6:Stream, title:str, out:str|None=None):
nt = st6[0].stats.npts
dt = st6[0].stats.delta
t = np.arange(nt)*dt

MAX, MIN = 0, 9e30
for tr in st6:
d = np.abs(tr.data)
if MAX < np.max(d):
MAX = np.max(d)
if MIN > np.min(d):
MIN = np.min(d)

travtP = stsyn[0].stats.sac['t0']
travtS = stsyn[0].stats.sac['t1']

fig, axs = plt.subplots(6, 1, figsize=(10, 7), gridspec_kw=dict(hspace=0.0), sharex=True)
for i in range(6):
tr = st6[i]
ax = axs[i]

ax.plot(t, tr.data, c='k', lw=0.5, label=tr.stats.channel)
ax.legend(loc='upper left')

ylims = ax.get_ylim()
if np.max(np.abs(np.array(ylims)))/MAX < 1e-5:
ylims = [-1, 1]
ax.set_ylim(ylims)

# 绘制到时
ax.vlines(travtP, *ylims, colors='b')
ax.text(travtP, ylims[1], "P", ha='left', va='top', color='b')
ax.vlines(travtS, *ylims, colors='r')
ax.text(travtS, ylims[1], "S", ha='left', va='top', color='r')

ax.set_xlim([t[0], t[-1]])
ax.set_ylim(np.array(ylims)*1.2)

axs[0].set_title(title)

if out is not None:
fig.savefig(out, dpi=100)



# BEGIN STRAIN
st_strain = pygrt.utils.compute_strain(stsyn)
print(st_strain)
# 6 Trace(s) in Stream:
# .SYN..ZZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..ZN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..ZE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..NN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..NE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..EE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# END STRAIN
plot6(st_strain, "Strain", "strain.png")

# BEGIN STRESS
st_stress = pygrt.utils.compute_stress(stsyn)
print(st_stress)
# 6 Trace(s) in Stream:
# .SYN..ZZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..ZN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..ZE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..NN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..NE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..EE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# END STRESS
plot6(st_stress, "Stress", "stress.png")
28 changes: 28 additions & 0 deletions docs/source/Tutorial/dynamic/run_upar/run.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
#!/bin/bash


# BEGIN GRN
# -e 表示计算空间导数
grt -Mmilrow -D2/0 -N500/0.02 -OGRN -R10 -e
# END GRN

# BEGIN SYN DC
# -e 表示计算空间导数
grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc -e
# END SYN DC

# BEGIN ZNE
# -N 表示输出ZNE分量
grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc_zne -e -N
# END ZNE


# BEGIN STRAIN
# 指定文件夹以及中间名，会在原文件夹内输出SAC格式的应变
grt.strain syn_dc_zne/out
# END STRAIN

# BEGIN STRESS
# 指定文件夹以及中间名，会在原文件夹内输出SAC格式的应变
grt.stress syn_dc_zne/out
# END STRESS
Binary file added docs/source/Tutorial/dynamic/run_upar/strain.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/Tutorial/dynamic/run_upar/stress.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading