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
3 changes: 2 additions & 1 deletion docs/source/Advanced/integ_converg/integ_converg.rst
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@
:start-after: BEGIN DGRN
:end-before: END DGRN

输出的核函数文件会在自定义路径下。
输出的核函数文件会在 :rst:dir:`GRN_grtstats/milrow_{depsrc}_{deprcv}/` 路径下
(与 C 一致,根目录名由 ``set_dynamic_grn_path`` 的输出目录决定)。


C和Python导出的核函数文件是一致的,底层调用的是相同的函数。文件名称格式为 ``K_{iw}_{freq}``,其中 ``{iw}`` 表示频率索引值, ``{freq}`` 表示对应频率(Hz)。文件为自定义的二进制文件, **强烈建议使用Python进行读取及后续处理**。这里还是给出两种读取方法。
Expand Down
3 changes: 2 additions & 1 deletion docs/source/Advanced/integ_converg/ptam.rst
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@
:start-after: BEGIN DEPSRC 0.0 DGRN
:end-before: END DEPSRC 0.0 DGRN

输出的核函数文件会在自定义路径下。
输出的核函数文件会在 :rst:dir:`GRN_grtstats/milrow_{depsrc}_{deprcv}/` 路径下
(与 C 一致)。

在 ``K_{iw}_{freq}`` 文件同级目录下,程序把 **PTAM过程中的核函数以及积分峰谷位置分为两个文件**
保存在 ``PTAM_{ir}_{dist}/`` 目录下( ``{ir}`` 为震中距索引, ``{dist}`` 为震中距),
Expand Down
34 changes: 23 additions & 11 deletions docs/source/Advanced/integ_converg/run/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,20 @@
# BEGIN DGRN
import numpy as np
import pygrt

modarr = np.loadtxt("milrow")
from pygrt.cli import format_float

depsrc = 2.0
deprcv = 0.0

pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)
pymod = pygrt.PyModel1D("milrow")
pymod.set_dynamic_grn_path("GRN")

# 通过statsfile参数自定义核函数文件的输出目录, statsidxs指定想输出的频率索引值
# statsidxs 指定频率索引,核函数写入 GRN_grtstats/{model}_{depsrc}_{deprcv}/
distarr = [5,8,10]
stgrnLst = pymod.compute_grn(
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv,
distarr=distarr, nt=500, dt=0.02,
statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
statsidxs=[50,100],
)
# END DGRN
# -------------------------------------------------------------------
Expand All @@ -24,7 +25,8 @@
# BEGIN read statsfile
# 可使用通配符简化输入,因为对应索引值下只会有一个文件
# 返回的是自定义类型的numpy数组
statsdata = pygrt.utils.read_statsfile(f"pygrtstats_{depsrc}_{deprcv}/K_0050_*")
statsdir = f"GRN_grtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}"
statsdata = pygrt.utils.read_statsfile(f"{statsdir}/K_0050_*")
print(statsdata.dtype)
# [('k', '<f8'), ('EX_q', '<c16'), ('EX_w', '<c16'), ('VF_q', '<c16'), ('VF_w', '<c16'), ('HF_q', '<c16'), ('HF_w', '<c16'), ('HF_v', '<c16'), ('DD_q', '<c16'), ('DD_w', '<c16'), ('DS_q', '<c16'), ('DS_w', '<c16'), ('DS_v', '<c16'), ('SS_q', '<c16'), ('SS_w', '<c16'), ('SS_v', '<c16')]
# END read statsfile
Expand Down Expand Up @@ -59,14 +61,17 @@
# BEGIN DEPSRC 0.0 DGRN
depsrc = 0.0
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)
pymod = pygrt.PyModel1D("milrow")
pymod.set_dynamic_grn_path("GRN")

stgrnLst = pymod.compute_grn(
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv,
distarr=distarr, nt=500, dt=0.02, converg_method='none',
statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
statsidxs=[50,100],
)

statsdata = pygrt.utils.read_statsfile(f"pygrtstats_{depsrc}_{deprcv}/K_0050_*")
statsdir = f"GRN_grtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}"
statsdata = pygrt.utils.read_statsfile(f"{statsdir}/K_0050_*")

dist=10
srctype="SS"
Expand All @@ -76,3 +81,10 @@
# END DEPSRC 0.0 DGRN
# -------------------------------------------------------------------

# 删除中间计算结果,仅保留成图
import shutil
from pathlib import Path
for name in ["GRN", "GRN_grtstats"]:
p = Path(name)
if p.is_dir():
shutil.rmtree(p, ignore_errors=True)
49 changes: 34 additions & 15 deletions docs/source/Advanced/integ_converg/run_ptam/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,18 +2,19 @@
# BEGIN DEPSRC 0.0 DGRN
import numpy as np
import pygrt

modarr = np.loadtxt("milrow")
from pygrt.cli import format_float

depsrc = 0.0
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)
pymod = pygrt.PyModel1D("milrow")
pymod.set_dynamic_grn_path("GRN")

distarr = [5,8,10]
# 设置 converg_method='PTAM' 进行收敛
stgrnLst = pymod.compute_grn(
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv,
distarr=distarr, nt=500, dt=0.02, converg_method='PTAM', k0=2, ampk=1.2, use_kmax_ref=True,
statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
statsidxs=[50,100],
)
# END DEPSRC 0.0 DGRN
# -------------------------------------------------------------------
Expand All @@ -22,7 +23,9 @@
# -------------------------------------------------------------------
# BEGIN plot ptam
ir = 2
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"pygrtstats_{depsrc}_{deprcv}/PTAM_{ir:04d}_*/PTAM_0050_*")
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(
f"GRN_grtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}/PTAM_{ir:04d}_*/PTAM_0050_*"
)

srctype="SS"
ptype="0"
Expand All @@ -36,20 +39,26 @@
# -------------------------------------------------------------------
# BEGIN SGRN
import numpy as np
import pygrt

modarr = np.loadtxt("milrow")
import pygrt
from pygrt.cli import format_float

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

norths = np.array([2.0])
easts = np.array([2.0])
static_grn = pymod.compute_static_grn(norths, easts, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}", k0=3, use_kmax_ref=True)
pymod = pygrt.PyModel1D("milrow")
pymod.set_static_grn_path("stgrn.nc")

norths = [2.0, 2.0, 1.0]
easts = [2.0, 2.0, 1.0]
pymod.compute_static_grn(
depsrc=depsrc, deprcv=deprcv,
norths=norths, easts=easts,
converg_method='PTAM', stats=True, k0=3, use_kmax_ref=True,
)

ir = 0
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"static_pygrtstats_{depsrc}_{deprcv}/PTAM_{ir:04d}_*/PTAM")
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(
f"stgrtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}/PTAM_{ir:04d}_*/PTAM"
)

srctype="SS"
ptype="0"
Expand All @@ -65,3 +74,13 @@

# END SGRN
# -------------------------------------------------------------------

# 删除中间计算结果,仅保留成图
import shutil
from pathlib import Path
for name in ["GRN", "GRN_grtstats", "stgrn.nc", "stgrtstats"]:
p = Path(name)
if p.is_dir():
shutil.rmtree(p, ignore_errors=True)
elif p.is_file():
p.unlink(missing_ok=True)
5 changes: 2 additions & 3 deletions docs/source/Advanced/integ_converg/run_ptam/run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,7 @@ rm -rf GRN* syn* ptam* pygrtstats* static* stgrtstats *.nc *.svg
# -------------------------------------------------------------------
# BEGIN DEPSRC 0.0 DGRN
# 使用 -Cp 指定使用 PTAM 进行收敛
grt greenfn -Mmilrow -D0/0 -N500/0.02 -OGRN -R5,8,10 -Cp -S50,100
grt ker2asc GRN_grtstats/milrow_0_0/K_0050_5.00000e+00 > stats
grt greenfn -Mmilrow -D0/0 -N500/0.02 -OGRN -R5,8,10 -Cp -K+k2+f+s1.2 -S50,100
# 绘制图像部分见Python
# END DEPSRC 0.0 DGRN
# -------------------------------------------------------------------
Expand All @@ -26,7 +25,7 @@ echo "..." >> ptam_stats_head
# -------------------------------------------------------------------
# BEGIN SGRN
# -S 表示输出核函数文件
grt static greenfn -Mmilrow -D0.05/0 -X2/2/1 -Y2/2/1 -Cp -S -Ostgrn.nc
grt static greenfn -Mmilrow -D0.05/0 -X2/2/1 -Y2/2/1 -Cp -K+k3+f -S -Ostgrn.nc

# grt.ker2asc 也可以读取静态解输出的核函数文件,格式一致
grt ker2asc stgrtstats/milrow_0.05_0/K > static_stats
Expand Down
31 changes: 25 additions & 6 deletions docs/source/Advanced/k_integ/drift/run/run.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,22 @@
# BEGIN 1
import numpy as np
import pygrt
from obspy import read
from pygrt.cli import format_float

modarr = np.loadtxt("milrow")

pymod = pygrt.PyModel1D(modarr, depsrc=10, deprcv=0.0)
pymod = pygrt.PyModel1D("milrow")
pymod.set_dynamic_grn_path("GRN")

depsrc = 10.0
deprcv = 0.0
nt = 500
dt = 10

st_grn = pymod.compute_grn(5000, nt=nt, dt=dt, keepAllFreq=True, statsfile="pygrtstats")[0]
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv, distarr=[5000],
nt=nt, dt=dt, keepAllFreq=True, statsidxs=[0, 1, 2, 3, 4, 5],
)
st_grn = read("GRN/*/*.sac")
# END 1

# 仅绘制一个分量做示例
Expand Down Expand Up @@ -85,7 +92,8 @@ def plot_freqs(tr:Trace):
# =================================================================
# 读入核函数
import glob
paths = glob.glob("pygrtstats/K_000[0-5]_*")
statsdir = f"GRN_grtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}"
paths = glob.glob(f"{statsdir}/K_000[0-5]_*")
paths.sort()
print(paths)

Expand All @@ -110,7 +118,10 @@ def plot_freqs(tr:Trace):

# =================================================================
# 跳过频段,重新计算
st_grn3 = pymod.compute_grn(5000, nt=nt, dt=dt)[0]
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv, distarr=[5000], nt=nt, dt=dt,
)
st_grn3 = read("GRN/*/*.sac")

srctypes = ['EX', 'VF', 'HF', 'DD', 'DS', 'SS']
chlst = ['Z', 'R', 'T']
Expand All @@ -137,3 +148,11 @@ def plot_all_waves(st_grn:Stream):
fig = plot_all_waves(st_grn3.copy())
fig.savefig("grn3.svg", bbox_inches='tight')
# =================================================================

# 删除中间计算结果,仅保留成图
import shutil
from pathlib import Path
for name in ["GRN", f"GRN_grtstats"]:
p = Path(name)
if p.is_dir():
shutil.rmtree(p, ignore_errors=True)
6 changes: 3 additions & 3 deletions docs/source/Advanced/k_integ/kmax.rst
Original file line number Diff line number Diff line change
Expand Up @@ -72,8 +72,8 @@
积分上限

:func:`compute_static_grn() <pygrt.pymod.PyModel1D.compute_static_grn>` 函数支持以下可选参数来控制波数积分,
参数与上面对应,具体说明详见API。
具体说明详见API。静态解无 ``ampk`` 参数

+ ``k0:float``
+ ``keps:float``
+ ``k0:float``
+ ``keps:float``
+ ``use_kmax_ref:bool``
26 changes: 17 additions & 9 deletions docs/source/Advanced/k_integ/safim/run/run.py
Original file line number Diff line number Diff line change
@@ -1,16 +1,24 @@
import numpy as np
import pygrt

modarr = np.loadtxt("milrow")
pymod = pygrt.PyModel1D("milrow")
pymod.set_dynamic_grn_path("GRN")

pymod = pygrt.PyModel1D(modarr, 5, 0)

st_grn = pymod.compute_grn(
distarr=[2500],
nt=2000,
dt=1,
Length=20,
pymod.compute_grn(
depsrc=5.0,
deprcv=0.0,
distarr=[2500],
nt=2000,
dt=1,
Length=20,
safilonTol=1e-2, # 自适应采样精度
filonCut=10,
delayT0=100,
)[0]
)

# 删除中间计算结果(成图由 plot.py 负责)
import shutil
from pathlib import Path
p = Path("GRN")
if p.is_dir():
shutil.rmtree(p, ignore_errors=True)
26 changes: 21 additions & 5 deletions docs/source/Advanced/kernel_old/run/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,22 @@
import matplotlib.pyplot as plt
from typing import Union
import pygrt
from pygrt.cli import format_float

modarr = np.loadtxt("mod1")
pymod = pygrt.PyModel1D("mod1")
pymod.set_dynamic_grn_path("KERNEL")

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

# 不指定statsidx,默认输出全部频率点的积分过程文件
# 不指定 statsidxs 索引时传空列表,输出全部频率点的积分过程文件
# vmin_ref 显式给定参考速度(用于定义波数积分上限),避免使用PTAM
# Length 给定波数积分间隔dk
_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, use_kmax_ref=True, converg_method='none', statsfile="pygrtstats")
pymod.compute_grn(
depsrc=depsrc, deprcv=deprcv, distarr=[1], nt=500, dt=0.02,
vmin_ref=0.1, Length=20, use_kmax_ref=True, converg_method='none',
statsidxs=[],
)
# END GRN
# -----------------------------------------------------------------

Expand All @@ -23,7 +30,8 @@

# 读取所有频率的核函数,并插值到vels
# 不指定ktypes,默认返回全部核函数,均以2D数组的形式保存,shape=(nfreqs, nvels)
kerDct = pygrt.utils.read_kernels_freqs("pygrtstats", vels)
statsdir = f"KERNEL_grtstats/mod1_{format_float(depsrc)}_{format_float(deprcv)}"
kerDct = pygrt.utils.read_kernels_freqs(statsdir, vels)
print(kerDct.keys())
# dict_keys(['_vels', '_freqs', 'EX_q', 'EX_w', 'VF_q', 'VF_w', 'HF_q', 'HF_w', 'HF_v', 'DD_q', 'DD_w', 'DS_q', 'DS_w', 'DS_v', 'SS_q', 'SS_w', 'SS_v'])
# END read
Expand Down Expand Up @@ -78,3 +86,11 @@ def plot_kernel(kerDct:dict, RorI:bool, out:Union[str,None]=None):
plot_kernel(kerDct, True, "real.svg")
# END plot
# -----------------------------------------------------------------

# 删除中间计算结果,仅保留成图
import shutil
from pathlib import Path
for name in ["KERNEL", "KERNEL_grtstats"]:
p = Path(name)
if p.is_dir():
shutil.rmtree(p, ignore_errors=True)
2 changes: 1 addition & 1 deletion docs/source/Advanced/kernel_old/run/run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ rm -rf GRN* pygrtstats* *.svg
# -S 后不指定索引表示输出所有频率点的核函数
# -Cn 禁用收敛算法
# -L20 定义波数积分间隔dk
grt greenfn -Mmod1 -D0.03/0 -N500/0.02 -OGRN -R1 -K+v0.1 -S -L20 -Cn
grt greenfn -Mmod1 -D0.03/0 -N500/0.02 -OGRN -R1 -K+v0.1+f -S -L20 -Cn
# END GRN
# -----------------------------------------------------------------

Expand Down
2 changes: 1 addition & 1 deletion docs/source/Gallery/ex01/ex01.rst
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
下载示例: :download:`ex01.tar.gz`

在 MILROW 模型下,震源深度 5 km,震源为剪切源,走向、倾角、滑动角分别为 77°、88°、99°,
计算震中距为 100 km ,方位角为 39.2 ° 的地面台站记录到的三分量理论记录。
计算震中距为 180 km ,方位角为 39.2 ° 的地面台站记录到的三分量理论记录。

此图也是 `PyGRT 代码主页 <https://github.com/Dengda98/PyGRT>`_
中显示的示例图。
Expand Down
Loading
Loading