Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
33 changes: 19 additions & 14 deletions docs/source/Gallery/ex13/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
from scipy.io import netcdf_file

def plot_disp(fLst):
fig, axs = plt.subplots(2, 3, figsize=(12, 6))
fig, axs = plt.subplots(3, 3, figsize=(12, 8))
for i, f in enumerate(fLst):
ax3 = axs[i]
chs = ['E', 'N', 'Z']
Expand All @@ -18,11 +18,12 @@ def plot_disp(fLst):

axs[0, 0].set_ylabel("Strike-slip", fontsize=12)
axs[1, 0].set_ylabel("Dip-slip", fontsize=12)
axs[2, 0].set_ylabel("Tensile", fontsize=12)
return fig, axs

def plot_stress(fLst):
fig, axs = plt.subplots(2, 2, figsize=(9,7), gridspec_kw=dict(hspace=0.2))
for i, f in enumerate([f1, f2]):
fig, axs = plt.subplots(2, 3, figsize=(12, 7), gridspec_kw=dict(hspace=0.2))
for i, f in enumerate(fLst):
ax = axs[0, i]
Y = f.variables['north'][:]
ax.plot(Y, f.variables['stress_ZZ'][:,0], 'bo-', ms=3, label='ZZ')
Expand All @@ -38,33 +39,37 @@ def plot_stress(fLst):

axs[0, 0].set_title("Strike-slip", fontsize=12)
axs[0, 1].set_title("Dip-slip", fontsize=12)
axs[0, 2].set_title("Tensile", fontsize=12)

axs[0, 0].set_ylabel("Normal", fontsize=12)
axs[1, 0].set_ylabel("Shear", fontsize=12)

axs[0, 1].legend(bbox_to_anchor=(1.3, 1))
axs[1, 1].legend(bbox_to_anchor=(1.3, 1))
axs[0, 2].legend(bbox_to_anchor=(1.3, 1))
axs[1, 2].legend(bbox_to_anchor=(1.3, 1))
return fig, axs


# ======================================================================
with netcdf_file("stsyn_1_stkslip.nc", mmap=False) as f1, \
netcdf_file("stsyn_1_dipslip.nc", mmap=False) as f2:
fig, axs = plot_disp([f1, f2])
netcdf_file("stsyn_1_dipslip.nc", mmap=False) as f2, \
netcdf_file("stsyn_1_tensile.nc", mmap=False) as f3:
fig, axs = plot_disp([f1, f2, f3])
fig.savefig("disp1.svg", bbox_inches='tight')
fig, axs = plot_stress([f1, f2])
fig, axs = plot_stress([f1, f2, f3])
fig.savefig("stress1.svg", bbox_inches='tight')

with netcdf_file("stsyn_2_stkslip.nc", mmap=False) as f1, \
netcdf_file("stsyn_2_dipslip.nc", mmap=False) as f2:
fig, axs = plot_disp([f1, f2])
netcdf_file("stsyn_2_dipslip.nc", mmap=False) as f2, \
netcdf_file("stsyn_2_tensile.nc", mmap=False) as f3:
fig, axs = plot_disp([f1, f2, f3])
fig.savefig("disp2.svg", bbox_inches='tight')
fig, axs = plot_stress([f1, f2])
fig, axs = plot_stress([f1, f2, f3])
fig.savefig("stress2.svg", bbox_inches='tight')

with netcdf_file("stsyn_3_stkslip.nc", mmap=False) as f1, \
netcdf_file("stsyn_3_dipslip.nc", mmap=False) as f2:
fig, axs = plot_disp([f1, f2])
netcdf_file("stsyn_3_dipslip.nc", mmap=False) as f2, \
netcdf_file("stsyn_3_tensile.nc", mmap=False) as f3:
fig, axs = plot_disp([f1, f2, f3])
fig.savefig("disp3.svg", bbox_inches='tight')
fig, axs = plot_stress([f1, f2])
fig, axs = plot_stress([f1, f2, f3])
fig.savefig("stress3.svg", bbox_inches='tight')
7 changes: 7 additions & 0 deletions docs/source/Gallery/ex13/run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@ grt static stress stsyn_1_stkslip.nc
grt static syn -M90/70/90 -Su1e6 -e -N -Gstgrn1.nc -Ostsyn_1_dipslip.nc
grt static stress stsyn_1_dipslip.nc

grt static syn -M90/70 -Su1e6 -e -N -Gstgrn1.nc -Ostsyn_1_tensile.nc
grt static stress stsyn_1_tensile.nc


grt static greenfn -M$mod -D10/14 -Y3/3/1 -X-10/10/0.5 -e -Ostgrn2.nc

Expand All @@ -22,6 +25,8 @@ grt static stress stsyn_2_stkslip.nc
grt static syn -M90/70/90 -Su1e6 -e -N -Gstgrn2.nc -Ostsyn_2_dipslip.nc
grt static stress stsyn_2_dipslip.nc

grt static syn -M90/70 -Su1e6 -e -N -Gstgrn2.nc -Ostsyn_2_tensile.nc
grt static stress stsyn_2_tensile.nc

mod="mod"
grt static greenfn -M$mod -D5/0 -Y3/3/1 -X-10/10/0.5 -e -Ostgrn3.nc
Expand All @@ -32,6 +37,8 @@ grt static stress stsyn_3_stkslip.nc
grt static syn -M90/70/90 -Su1e6 -e -N -Gstgrn3.nc -Ostsyn_3_dipslip.nc
grt static stress stsyn_3_dipslip.nc

grt static syn -M90/70 -Su1e6 -e -N -Gstgrn3.nc -Ostsyn_3_tensile.nc
grt static stress stsyn_3_tensile.nc


python plot.py
Expand Down
Empty file modified docs/source/Gallery/ex14/run.sh
100644 → 100755
Empty file.
17 changes: 17 additions & 0 deletions docs/source/Tutorial/dynamic/run/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,23 @@ def plot_int_dif(stsyn:Stream, stsyn_int:Stream, stsyn_dif:Stream, comp:str, out
# -----------------------------------------------------------------------------------


# -----------------------------------------------------------------------------------
# BEGIN SYN TS
# 接之前的代码
idx = 2
stgrn = stgrnLst[idx] # 选择格林函数

stsyn = pygrt.utils.gen_syn_from_gf_TS(stgrn, M0=1e24, strike=33, dip=50, az=30)
print(stsyn)
# 3 Trace(s) in Stream:
# .SYN..Z | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..R | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
# .SYN..T | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples
plot_syn(stsyn, "syn_ts.svg")
# END SYN TS
# -----------------------------------------------------------------------------------


# -----------------------------------------------------------------------------------
# BEGIN SYN MT
idx = 2
Expand Down
7 changes: 7 additions & 0 deletions docs/source/Tutorial/dynamic/run/run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,13 @@ grt syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc
# END SYN DC
# -----------------------------------------------------------------------------------

# -----------------------------------------------------------------------------------
# BEGIN SYN TS
# 合成结果在 syn_ts/ 目录下,以SAC格式保存。
grt syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50 -Osyn_ts
# END SYN TS
# -----------------------------------------------------------------------------------

# -----------------------------------------------------------------------------------
# BEGIN SYN MT
# 合成结果在 syn_mt/ 目录下,以SAC格式保存。
Expand Down
25 changes: 25 additions & 0 deletions docs/source/Tutorial/dynamic/syn.rst
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,31 @@ Python中合成动态位移的主函数为 :func:`gen_syn_from_gf_*() <pygrt.uti
:align: center


张裂源
~~~~~~~~~~~~~~
断层走向33°,倾角50°,标量矩 1e24 dyne·cm。

.. tabs::

.. group-tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN SYN TS
:end-before: END SYN TS

.. group-tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN SYN TS
:end-before: END SYN TS


.. figure:: run/syn_ts.svg
:align: center


矩张量源
~~~~~~~~~~~~~~
:math:`M_{xx}=0.1, M_{xy}=-0.2, M_{xz}=1.0, M_{yy}=0.3, M_{yz}=-0.5, M_{zz}=-2.0`,单位 1e24 dyne·cm, **其中X为北向,Y为东向,Z为垂直向下**。
Expand Down
18 changes: 18 additions & 0 deletions docs/source/Tutorial/static/run/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,24 @@ def plot_static(static_syn:dict, out:Union[str,None]=None):
# END SYN DC2
# ---------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------
# BEGIN SYN TS
static_syn = pygrt.utils.gen_syn_from_gf_TS(static_grn, M0=1e24, strike=33, dip=50, ZNE=True)
print(static_syn.keys())
# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'Z', 'N', 'E'])
plot_static(static_syn, "syn_ts.svg")
# END SYN TS
# ---------------------------------------------------------------------------------


# ---------------------------------------------------------------------------------
# BEGIN SYN TS2
static_syn = pygrt.utils.gen_syn_from_gf_TS(static_grn, M0=1e24, strike=33, dip=90, ZNE=True)
print(static_syn.keys())
# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'Z', 'N', 'E'])
plot_static(static_syn, "syn_ts2.svg")
# END SYN TS2
# ---------------------------------------------------------------------------------


# ---------------------------------------------------------------------------------
Expand Down
33 changes: 33 additions & 0 deletions docs/source/Tutorial/static/run/run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,39 @@ EOF
gmt colorbar -Bx+l"Z (cm)"
gmt end

# ---------------------------------------------------------------------------------
# BEGIN SYN TS
# 从网格文件中读取格林函数,再将合成结果写入新网格
grt static syn -S1e24 -M33/50 -N -Gstgrn.nc -Ostsyn_ts.nc
# END SYN TS
# ---------------------------------------------------------------------------------

# 计算张裂的矩张量表示,仅绘制其中的 DC+CLVD 分量
gmt begin syn_ts pdf
gmtplot_static stsyn_ts.nc -Si0.03c
gmt meca -Sz0.5c <<EOF
0 0 2 $(python tension2mt.py 33 50 stsyn_ts.nc) 24
EOF
gmt colorbar -Bx+l"Z (cm)"
gmt end

# ---------------------------------------------------------------------------------
# BEGIN SYN TS2
# 从网格文件中读取格林函数,再将合成结果写入新网格
grt static syn -S1e24 -M33/90 -N -Gstgrn.nc -Ostsyn_ts2.nc
# END SYN TS2
# ---------------------------------------------------------------------------------

# 计算张裂的矩张量表示,仅绘制其中的 DC+CLVD 分量
gmt begin syn_ts2 pdf
gmtplot_static stsyn_ts2.nc -Si0.03c
gmt meca -Sz0.5c <<EOF
0 0 2 $(python tension2mt.py 33 90 stsyn_ts2.nc) 24
EOF
gmt colorbar -Bx+l"Z (cm)"
gmt end


# ---------------------------------------------------------------------------------
# BEGIN SYN MT
# 从网格文件中读取格林函数,再将合成结果写入新网格
Expand Down
34 changes: 34 additions & 0 deletions docs/source/Tutorial/static/run/tension2mt.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
"""
Zhu Dengda

将张裂转为矩张量表示
"""
from math import *
import numpy as np
from scipy.io import netcdf_file
import sys

strike = radians(float(sys.argv[1]))
dip = radians(float(sys.argv[2]))
ncfile = sys.argv[3]

nvec = np.array([
- sin(dip) * sin(strike),
sin(dip) * cos(strike),
cos(dip),
])

vvec = nvec.copy()

# 从 nc 文件中读取 src_va, src_vb
with netcdf_file(ncfile, mmap=False) as f:
src_va = f._attributes['src_va']
src_vb = f._attributes['src_vb']

M = ((src_va/src_vb)**2 - 2) * np.eye(3) * np.sum(nvec * vvec) + (np.einsum('i,j', nvec, vvec) + np.einsum('i,j', vvec, nvec))

# 输出为 mrr mtt mff mrt mrf mtf
Mrtf = [
M[2,2], M[0,0], M[1,1], M[0,2], -M[1,2], -M[0,1]
]
print(*Mrtf, sep=' ', end='')
52 changes: 52 additions & 0 deletions docs/source/Tutorial/static/static_syn.rst
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,58 @@ Python中合成静态位移的主函数为 :func:`gen_syn_from_gf_*() <pygrt.uti
:align: center


张裂源
~~~~~~~~~~~~~~
断层走向33°,倾角50°,标量矩 1e24 dyne·cm。

.. tabs::

.. group-tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN SYN TS
:end-before: END SYN TS

.. group-tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN SYN TS
:end-before: END SYN TS


.. figure:: run/syn_ts.svg
:width: 500px
:align: center

沙滩球仅绘制了张裂源中的 DC+CLVD 分量


这里如果改变倾角为90°,就可以看到清晰的对称辐射花样。

.. tabs::

.. group-tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN SYN TS2
:end-before: END SYN TS2

.. group-tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN SYN TS2
:end-before: END SYN TS2


.. figure:: run/syn_ts2.svg
:width: 500px
:align: center

沙滩球仅绘制了张裂源中的 DC+CLVD 分量


矩张量源
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/radiation.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
X(GRT_SYN_EX, "Explosion") \
X(GRT_SYN_SF, "Single Force") \
X(GRT_SYN_DC, "Shear") \
X(GRT_SYN_TS, "Tensile") \
X(GRT_SYN_TS, "Tension") \
X(GRT_SYN_MT, "Moment Tensor") \


Expand All @@ -44,7 +44,7 @@ static const char *srcTypeFullName[] = {
* @param[in] par_theta 方向因子中是否对theta(az)求导
* @param[in] M0 放大系数,对于剪切源、爆炸源、张量震源,M0是标量地震矩;对于单力源,M0是放大系数
* @param[in] coef 放大系数,用于位移空间导数的计算
* @param[in] VpVs_ratio 震源层的 Vp/Vs 比值,用于张位错源
* @param[in] VpVs_ratio 震源层的 Vp/Vs 比值,用于张裂源
* @param[in] azrad 弧度制的方位角
* @param[in] mchn 震源机制参数,
* 对于单力源,mchn={fn, fe, fz},
Expand Down
2 changes: 1 addition & 1 deletion pygrt/C_extension/src/dynamic/grt_syn.c
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
};
break;

// 剪切震源, 张位错源
// 剪切震源, 张裂源
case 'M':
Ctrl->M.active = true;
{
Expand Down
2 changes: 1 addition & 1 deletion pygrt/C_extension/src/static/grt_static_syn.c
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
};
break;

// 剪切震源, 张位错源
// 剪切震源, 张裂源
case 'M':
Ctrl->M.active = true;
{
Expand Down
6 changes: 3 additions & 3 deletions pygrt/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
:date: 2024-07-24

该文件包含一些数据处理操作上的补充:
1、剪切源、张位错源、单力源、爆炸源、矩张量源 通过格林函数合成理论地震图的函数\n
1、剪切源、张裂源、单力源、爆炸源、矩张量源 通过格林函数合成理论地震图的函数\n
2、Stream类型的时域卷积、微分、积分 (基于numpy和scipy) \n
3、Stream类型写到本地sac文件,自定义名称 \n
4、读取波数积分和峰谷平均法过程文件 \n
Expand Down Expand Up @@ -96,7 +96,7 @@ def _gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:fl
srcName = ["EX", "VF", "HF", "DD", "DS", "SS"]
allchs = [tr.stats.channel for tr in st]

# 为张位错计算 Vp/Vs
# 为张裂计算 Vp/Vs
src_va = st[0].stats.sac['user6']
src_vb = st[0].stats.sac['user7']
VpVs_ratio = src_va / src_vb
Expand Down Expand Up @@ -168,7 +168,7 @@ def _gen_syn_from_static_gf(grn:dict, calc_upar:bool, compute_type:GRT_SYN_TYPE,
srcName = ["EX", "VF", "HF", "DD", "DS", "SS"]
allchs = list(grn.keys())

# 为张位错计算 Vp/Vs
# 为张裂计算 Vp/Vs
src_va = grn['_src_va']
src_vb = grn['_src_vb']
VpVs_ratio = src_va / src_vb
Expand Down
Loading