pygrt.pymod

file:: pymod.py
author:: Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
date:: 2024-07-24

该文件包括 Python 端使用的基于文件的模型 PyModel1D

基类：object

File-based 1D layered model for GRT calculations.

Typical workflow:

Create PyModel1D with the Green’s function path(s) actually needed (grn and/or stgrn) and optional modelpath.
Compute modal dispersion and Green’s functions with eigenv(), eigenfn(), modsum(), greenfn() or static_greenfn() (the required paths depend on the method).
Synthesize waveforms or static fields with syn() or static_syn() (only the corresponding GF path is required).

Create a file-based 1D layered model handle.

The model file is a plain text table. Each row is one layer in the form thickness(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) [Qp Qs]. A zero thickness marks a half-space bottom layer.

All arguments must be passed by keyword. grn / stgrn / modelpath are optional at construction; methods that need them raise if missing.

参数:

grn – Root directory for dynamic Green’s functions.
stgrn – NetCDF file path for static Green’s functions.
modelpath – Path to the layered model file. Required by eigenv(), greenfn() and static_greenfn(). Travel-time calculations also need it, unless a per-call modelpath is passed to travt().
topbound – Top boundary condition. One of free, rigid and halfspace.
botbound – Bottom boundary condition. One of free, rigid and halfspace.

travt(*, depsrc: float, deprcv: float, dists: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], modelpath: str | PathLike | None = None)[源代码]

Compute first-arrival P- and S-wave travel times.

Calls the C routine grt_compute_travt1d_from_file, which reads the layered model and evaluates travel times at the given source/receiver depths and epicentral distances. All arguments must be passed by keyword.

参数:

depsrc – Source depth in km.
deprcv – Receiver depth in km.
dists – Epicentral distance(s) in km. A scalar or a strictly ascending sequence of distances.
modelpath – Model file for this call only. If omitted, uses self.modelpath. If both are set and differ, a warning is issued and self.modelpath is not changed. Required when self.modelpath is unset.

返回:

(travtP, travtS) in s. For a scalar distance both are floats; for multiple distances both are NumPy arrays of shape (n,).

compute_travt1d(*args, **kwargs)[源代码]: Legacy interface renamed to travt(); calling it raises an error.

Compute surface-wave phase-velocity dispersion with grt eigenv.

The ordinary workflow requires modelpath in the constructor and writes a NetCDF dispersion file to phase_path. Set secular_freq to use the C module’s debug mode, which prints the secular function instead and does not require freqs / periods or phase_path. All arguments must be passed by keyword.

参数:

wtype – Surface-wave type, "R" for Rayleigh or "L" for Love.
freqs – One frequency, or (f1, f2, df) in Hz. Mutually exclusive with periods.
periods – One period, or (T1, T2, dT) in s. Mutually exclusive with freqs.
phase_path – Output NetCDF dispersion path. It is not used in secular-function debug mode.
max_order – Maximum mode order. None leaves -N unspecified and computes the fundamental mode.
all_modes – If true, pass bare -N and compute as many modes as possible. Mutually exclusive with max_order.
secular_freq – Debug frequency in Hz. When set, pass -X instead of -F and -C.
cmin – Lower phase-velocity bound in km/s for debug secular-function output.
cmax – Upper phase-velocity bound in km/s for debug secular-function output.
iref – Layer index for debug secular-function output.
ctrl_kw – Optional dict of search-control parameters for -T. Allowed keys are tol, cgap, thrd, vgap and dc. tol, cgap, thrd and vgap are decimal exponents of the CLI thresholds (tol=3 means \(10^{-3}\)), while dc is a uniform search interval in km/s.
nthreads – Number of OpenMP threads. None uses all cores (CLI -P).
print_log – If false, suppress command output and capture it on failure. Warnings and diagnostics from grt are still printed.

返回:

None. The dispersion file is written to phase_path in ordinary mode.

Compute modal eigenfunctions and derived quantities with grt eigenfn.

The model path is read from the NetCDF dispersion file generated by eigenv(). Set one or more output paths to request eigenfunctions, group velocity, energy integrals, or phase/group-velocity sensitivity kernels. All arguments must be passed by keyword.

参数:

phase_path – Input NetCDF phase-velocity dispersion path.
freqs – Optional one-, two- or three-value frequency selection in Hz. Mutually exclusive with periods.
periods – Optional one-, two- or three-value period selection in s. Mutually exclusive with freqs.
modes – Optional one-, two- or three-value mode selection (n1[, n2[, dn]]).
all_modes – If true, pass bare -N. Mutually exclusive with modes.
eigenfn_path – Output NetCDF eigenfunction path.
depths – Optional one- or three-value depth selection for eigenfn_path in km.
group_path – Output NetCDF group-velocity dispersion path.
csens_path – Output phase-velocity sensitivity path.
usens_path – Output group-velocity sensitivity path.
egy_path – Output energy-integral path.
sensitivity_dz – Depth interval in km for sensitivity and energy-integral calculations.
print_log – If false, suppress command output and capture it on failure.

返回:

None. Requested NetCDF files are written to their output paths.

Compute surface-wave Green’s functions with grt modsum.

phase_path must be a dispersion file generated by eigenv() with an equally spaced frequency grid. Source and receiver depths may be scalars or strictly ascending sequences; all depth combinations are computed when sequences are supplied. Results are written to output_path or, when it is omitted, to the grn root configured in the constructor. All arguments must be passed by keyword.

参数:

phase_path – Input NetCDF phase-velocity dispersion path.
depsrc – Source depth or strictly ascending source-depth sequence in km.
deprcv – Receiver depth or strictly ascending receiver-depth sequence in km.
dists – Epicentral distance or strictly ascending distance sequence in km.
output_path – Output Green’s-function root. If omitted, use self.grn.
freqband – Optional (f1, f2) frequency range in Hz.
modes – Optional one-, two- or three-value mode selection (n1[, n2[, dn]]).
all_modes – If true, pass bare -N. Mutually exclusive with modes.
upsampling_n – Frequency-domain zero-padding factor.
delayT0 – Time delay at zero distance in s. Supplying this, delayV0 or ref_first_p enables -E.
delayV0 – Reference velocity for the time delay in km/s. Used only when ref_first_p is false.
ref_first_p – Whether to use the first P-wave arrival as the reference for the time delay (CLI -Ep).
gf_source – Source types to calculate from EX, VF, HF and DC. None means all types.
calc_upar – If true, also compute spatial displacement derivatives.
nthreads – Number of OpenMP threads. None uses all cores (CLI -P).
print_log – If false, suppress command output and capture it on failure. Warnings and diagnostics from grt are still printed.

返回:

None. Green’s-function SAC files are written below the output root.

greenfn(*, depsrc: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], deprcv: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], dists: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], nt: int, dt: float, upsampling_n: int = 1, freqband: Sequence[float] = (-1.0, -1.0), zeta: float = 0.8, keepAllFreq: bool = False, vmin_ref: float = 0.0, keps: float = -1.0, ampk: float = 2.0, k0: float = 50.0, use_kmax_ref: bool = False, Length: float = 0.0, filonLength: float = 0.0, safilonTol: float = 0.0, filonCut: float = 0.0, converg_method: str = 'AUTO', delayT0: float = 0.0, delayV0: float = 0.0, ref_first_p: bool = False, skipImagComps: bool = False, calc_upar: bool = False, gf_source: Iterable[str] | None = None, statsidxs: Sequence[int] | None = None, nthreads: int | None = None, print_log: bool = True)[源代码]

Compute dynamic Green’s functions with the grt greenfn command.

Requires grn and modelpath in the constructor. Results are written as SAC files under {grn}/{model}_{depsrc}_{deprcv}_{distance}. The model file is also copied into the grn root directory by the C module. All arguments must be passed by keyword.

参数:

depsrc – Source depth or strictly ascending source-depth sequence in km.
deprcv – Receiver depth or strictly ascending receiver-depth sequence in km.
dists – Array of strictly ascending epicentral distances in km, or a single distance.
nt – Number of time points. With the help of SciPy, nt no longer needs to be a power of 2.
dt – Time interval in s.
upsampling_n – Upsampling factor applied after inverse FFT.
freqband – Frequency range (f1, f2) in Hz. Negative values mean that the corresponding bound is determined automatically.
zeta – Coefficient defining the imaginary angular frequency, \(\tilde{\omega} = \omega - j w_I\), where \(w_I = \zeta\pi/T\) and \(T=nt\,dt\).
keepAllFreq – Whether to calculate all frequency points, regardless of how low the frequency is.
vmin_ref – Minimum reference velocity in km/s. 0.0 means the minimum model velocity, limited to 0.1 km/s.
keps – Automatic convergence condition. See Yao and Harkrider (1983) for more details. A negative value disables this condition.
ampk – Amplification factor in the reference maximum wavenumber.
k0 – Coefficient in the reference maximum wavenumber, \(k_{\mathrm{max,ref}} = \sqrt{(k_0\pi/h_s)^2 + (ampk\,\omega/v_{\mathrm{min,ref}})^2}\). Here \(h_s=\max(|depsrc-deprcv|,0.1)\).
use_kmax_ref – Whether to use the reference maximum wavenumber directly, without amplitude searching.
Length – Integration step \(dk=2\pi/(Lr_{\max})\). 0.0 means the value is selected automatically.
filonLength – Integration step of fixed-interval Filon integration at large distances, but not at zero distance. 0.0 disables this method. Do not set together with safilonTol.
safilonTol – Precision of self-adaptive Filon integration at large distances, but not at zero distance. 0.0 disables this method. Do not set together with filonLength.
filonCut – Splitting point between DWM and Filon integration, \(k^*=\mathrm{filonCut}/r_{\max}\).
converg_method – Explicit convergence method. One of AUTO, NONE, DCM and PTAM.
delayT0 – Time delay at zero distance in s.
delayV0 – Reference velocity for the time delay in km/s. Used only when ref_first_p is false.
ref_first_p – Whether to use the first P-wave arrival as the reference for the time delay (CLI -Ep).
skipImagComps – Whether to skip the amplitude compensation from the imaginary frequency.
calc_upar – Whether to calculate spatial derivatives of displacement. Required later if strain, stress or rotation will be computed.
gf_source – Source types to calculate. Choose from EX, VF, HF and DC. None means all available source types.
statsidxs – Frequency indexes for optional statistics output. None means no statistics files. An empty list means all frequency indexes (CLI bare -S).
nthreads – Number of OpenMP threads. None uses all cores (CLI -P).
print_log – Whether to print calculation logs. Warnings and diagnostics from grt are always printed.

返回:

None. Results are written to disk.

compute_grn(*args, **kwargs)[源代码]: Legacy interface renamed to greenfn(); calling it raises an error.

static_greenfn(*, depsrc: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], deprcv: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]], dists: float | Sequence[float] | ndarray[Any, dtype[floating[Any]]] | None = None, norths: Sequence[float] | None = None, easts: Sequence[float] | None = None, keps: float = -1.0, k0: float = 50.0, use_kmax_ref: bool = False, Length: float = 15.0, filonLength: float = 0.0, safilonTol: float = 0.0, filonCut: float = 0.0, converg_method: str = 'AUTO', calc_upar: bool = False, stats: bool = False)[源代码]

Compute static Green’s functions with the grt static_greenfn command.

Requires stgrn and modelpath in the constructor. Results are written to the configured NetCDF file (4D STGRNLIB layout [depsrc][deprcv][north][east]) and currently overwrite any existing content. All arguments must be passed by keyword.

depsrc / deprcv may be a scalar or a 1-D sequence:

Single depth pair: CLI -Ddepsrc/deprcv.
Multiple depths: CLI -Ds... / -Dr... (comma-separated list).

For an ordinary Green’s-function library, dists is the usual choice: it computes a one-dimensional distance list and stores it as north=0 and east=dists. This is also the most convenient format for static_syn(), which synthesizes results at surrounding distance samples and combines those results with distance-based weights. Receiver locations can be specified in either of two ways:

dists, a scalar or sequence of epicentral distances in km. This is equivalent to placing receivers along the east axis with north = 0 and is the recommended way to build a reusable library.
norths and easts, each a three-value sequence (start, stop, step) in km, mapped to CLI -X / -Y. This form is mainly useful for accuracy tests or when a later synthesis should reuse the same grid without distance-based combination of multiple synthesized results.

参数:

depsrc – Source depth(s) in km. Multiple values must be strictly ascending.
deprcv – Receiver depth(s) in km. Multiple values must be strictly ascending.
dists – Epicentral distance(s) in km. A scalar or strictly ascending sequence, equivalent to receivers with north = 0 and east = dists. Mutually exclusive with norths / easts.
norths – Three values defining the north-coordinate option -Xstart/stop/step in km.
easts – Three values defining the east-coordinate option -Ystart/stop/step in km.
keps – Automatic convergence condition. See Yao and Harkrider (1983) for more details. A negative value disables this condition.
k0 – Coefficient in the reference maximum wavenumber, where \(h_s=\max(|depsrc-deprcv|,0.1)\). The actual maximum wavenumber is searched based on the kernel amplitude.
use_kmax_ref – Whether to use the reference maximum wavenumber directly, without amplitude searching.
Length – Integration step \(dk=2\pi/(Lr_{\max})\). The default is 15.
filonLength – Step parameter for fixed-interval Filon integration at large distances. 0.0 disables this method. Do not set together with safilonTol.
safilonTol – Tolerance for self-adaptive Filon integration at large distances. 0.0 disables this method. Do not set together with filonLength.
filonCut – Splitting point between DWM and Filon integration.
converg_method – Explicit convergence method. One of AUTO, NONE, DCM and PTAM.
calc_upar – Whether to calculate spatial derivatives of displacement. Required later if strain, stress or rotation will be computed.
stats – Whether to write integration statistics. Only available for a single source/receiver depth; ignored with a warning for multi-depth runs.

返回:

None. Results are written to the configured NetCDF file.

compute_static_grn(*args, **kwargs)[源代码]: Legacy interface renamed to static_greenfn(); calling it raises an error.

syn(*, depsrc: float | None = None, deprcv: float | None = None, dist: float | None = None, azimuth: float, output_path: str | PathLike, scale: float, scale_with_mu: bool = False, strike: float | None = None, dip: float | None = None, rake: float | None = None, force: Sequence[float] | None = None, moment_tensor: Sequence[float] | None = None, time_function: str | None = None, integrate_order: int | None = None, differentiate_order: int | None = None, zne: bool = False, calc_upar: bool = False, return_result: bool = False)[源代码]

Synthesize dynamic three-component displacement with grt syn.

Results are written as SAC files under output_path. By default the synthetics are impulse-like displacements in cm with ZRT components:

Z - vertical upward
R - radial outward
T - clockwise 90° from R

Call greenfn() first (or point grn at an existing GF root or subdirectory). When grn is a root, a selector is required for a dimension with multiple values. For a singleton dimension it may be omitted or explicitly set, but an explicit value must match. When grn is a subdirectory, all three selectors must be omitted. No interpolation is performed. All arguments must be passed by keyword.

The source type is inferred from the source-specific parameters. Leave strike, dip, rake, force and moment_tensor unset for an explosion (EX). Supplying force selects a single force (SF); supplying moment_tensor selects a moment tensor (MT); supplying strike and dip selects a tensile crack (TS), or a double-couple (DC) when rake is also supplied. Only one source parameter group may be used at a time.

参数:

depsrc – Source depth in km when grn is a GF root with multiple source depths. It may be omitted when the root has one source depth, but an explicit value must match the root library. It must be omitted when grn is a subdirectory.
deprcv – Receiver depth in km when grn is a GF root with multiple receiver depths. It may be omitted when the root has one receiver depth, but an explicit value must match the root library. It must be omitted when grn is a subdirectory.
dist – Epicentral distance in km when grn is a GF root with multiple distances. It may be omitted when the root has one distance, but an explicit value must match the root library. It must be omitted when grn is a subdirectory.
azimuth – Azimuth from source to receiver in deg. North is 0°, clockwise positive.
output_path – Output directory for SAC files {output_path}/{ch}.sac.
scale – Source scaling factor. For explosion, double-couple, tensile-crack and moment-tensor sources, this is the scalar seismic moment in dyne·cm. For a single force, the unit is dyne. If scale_with_mu is true, scale is treated as area × slip in cm³ and multiplied by the source-layer shear modulus \(\mu\).
scale_with_mu – If true, multiply scale by the source-layer shear modulus \(\mu\) (CLI -Su).
strike – Fault strike in deg, in [0, 360]. North is 0°, clockwise positive. Set together with dip for TS or with dip and rake for DC.
dip – Fault dip in deg, in [0, 90]. Required for TS and DC when the corresponding source geometry is selected.
rake – Slip rake in deg, in [-180, 180], counterclockwise positive on the fault plane. Supplying it with strike and dip selects DC; omit it for TS.
force – Single-force coefficients (fN, fE, fZ) for the SF source. Directions are north, east and downward. Each coefficient is multiplied by scale.
moment_tensor – Six independent moment-tensor coefficients (Mxx, Mxy, Mxz, Myy, Myz, Mzz) for the MT source. Subscripts x/y/z denote north/east/down.
time_function –
Time-function string passed to CLI -D. Supported forms include: * p/t0 - parabola lasting t0 s * t/t1/t2/t3 - trapezoid with rise,

plateau and fall cutoffs in s
- r/f0 - Ricker wavelet with dominant frequency f0 in Hz
- 0/file - custom file with exactly one amplitude column
All time functions use area normalization, except the Ricker wavelet, which has a peak amplitude of 1. A custom sequence should sum to 1; it is not normalized by the program, which only issues a warning otherwise. Omit this argument for an impulse response.
integrate_order – Number of time integrations. For example, 1 yields step-like displacement.
differentiate_order – Number of time differentiations. For example, 1 yields velocity.
zne – If true, output ZNE instead of ZRT components.
calc_upar – If true, also synthesize spatial derivatives of displacement. Derivative channel names are prefixed with z, r or t. Set this when strain, stress or rotation will be computed later.
return_result – If true, read the generated SAC files into an obspy.Stream.

返回:

An ObsPy stream when return_result is true; otherwise None.

compute_syn(*args, **kwargs)[源代码]: Legacy interface renamed to syn(); calling it raises an error.

Synthesize static three-component displacement with grt static_syn.

Results are written to the NetCDF file output_path. Requires stgrn (and typically a prior static_greenfn()). All arguments must be passed by keyword.

Receivers default to the library north/east grid. Optionally redefine them with norths/easts (uniform deprcv when the library has multiple receiver depths), or with recv_points for an ASCII file of arbitrary north east depth points (CLI -Q). Each row may append strike dip rake in degrees; these angles are saved in the output but are not used in synthesis. recv_points is mutually exclusive with norths/easts and deprcv. If the library was built with dists / -R, the default grid is a 1-D line (north = 0, east = R); set norths/easts or recv_points to obtain a 2-D field. A Coulomb-format finite receiver-fault file can be supplied through rcv_fault (CLI -R); without rcv_fault_size the subdivision size defaults to the smallest positive sampling interval among epicentral distance, source depth and receiver depth in the library. With rcv_fault_size, each fault is subdivided with (dL, dW) along strike / dip.

Point-source type is inferred from the source-specific parameters: leaving strike, dip, rake, force and moment_tensor unset selects an explosion (EX); force selects a single force (SF); moment_tensor selects a moment tensor (MT); strike and dip select a tensile crack (TS), or a double-couple (DC) when rake is also supplied. Only one source parameter group may be used at a time. Finite faults use src_fault (Coulomb-format file, CLI -C) instead. Its Kode column selects rectangular shear/ tensile sources or point shear/expansion sources. An exact rake token in the seventh header column selects Kode 100 rake/net-slip rows; the filename suffix is not used to select the format. That path requires a multi-source-depth library and rejects point-source options.

For each target receiver, the C module first synthesizes results at the surrounding epicentral-distance samples and combines those synthesized results with weights based on the target distance. When a requested source or receiver depth lies between samples, it performs the same process for each surrounding depth combination and then combines those synthesized results with depth-based weights. Thus, interpolation is applied to synthesized results rather than directly to Green’s-function arrays. If the library was generated with an explicit -X/-Y grid and the same grid is reused, the corresponding synthesized results can be used directly.

参数:

depsrc – Point-source depth in km (CLI -Ds). Required when the library has multiple source depths; optional when it has one, but an explicit value must match the library. Forbidden when src_fault is set.
deprcv – Receiver depth in km for grid receivers (CLI -Dr). Required when the library has multiple receiver depths and recv_points is not used; optional when it has one, but an explicit value must match the library. Do not set it when using recv_points.
norths – Optional new north grid as three values (start, stop, step) in km. Must be set together with easts. Mutually exclusive with recv_points.
easts – Optional new east grid as three values (start, stop, step) in km. Must be set together with norths. Mutually exclusive with recv_points.
recv_points – ASCII file of arbitrary receivers (north east depth in km, optionally followed by strike dip rake in degrees; # comments). All data rows must use the same 3- or 6-column format. Mutually exclusive with norths/easts and deprcv.
rcv_fault – Coulomb-format finite receiver-fault file with 11 data columns; an exact rake token in the seventh header column selects Kode 100 rake/net-slip interpretation (CLI -R). Without rcv_fault_size, the library sampling intervals determine the default subdivision size. With that argument, each fault contributes multiple subfault centers. Mutually exclusive with recv_points, norths/easts and deprcv.
rcv_fault_size – Optional positive (dL, dW) in km for receiver-fault subdivision along strike / dip; if omitted, use the smallest positive interval among epicentral distance, source depth and receiver depth in the library.
output_path – Output NetCDF file path.
scale – Point-source scaling factor. For explosion, double-couple, tensile-crack and moment-tensor sources, this is the scalar seismic moment in dyne·cm. For a single force, the unit is dyne. If scale_with_mu is true, scale is treated as area × slip in cm³ and multiplied by the source-layer shear modulus \(\mu\). Required for point sources; ignored for src_fault.
scale_with_mu – If true, multiply scale by the source-layer shear modulus \(\mu\) (CLI -Su).
strike – Fault strike in deg, in [0, 360]. North is 0°, clockwise positive. Set together with dip for TS or with dip and rake for DC.
dip – Fault dip in deg, in [0, 90]. Required for TS and DC when the corresponding source geometry is selected.
rake – Slip rake in deg, in [-180, 180], counterclockwise positive on the fault plane. Supplying it with strike and dip selects DC; omit it for TS.
force – Single-force coefficients (fN, fE, fZ) for the SF source. Directions are north, east and downward. Each coefficient is multiplied by scale.
moment_tensor – Six independent moment-tensor coefficients (Mxx, Mxy, Mxz, Myy, Myz, Mzz) for the MT source. Subscripts x/y/z denote north/east/down.
src_fault – Coulomb-format finite-fault file with 11 data columns (Kode 100/200/300/400/500). An exact rake token in the seventh header column selects Kode 100 rake/net-slip interpretation; the filename suffix is not used. Mutually exclusive with point-source options; point-source arguments cause ValueError.
src_fault_size – Optional (dL, dW) in km for finite-fault subdivision along strike / dip. If omitted, the C code uses the smallest positive interval among epicentral distance, source depth and receiver depth in the library.
zne – If true, output ZNE instead of ZRT components.
calc_upar – If true, also synthesize spatial derivatives of displacement. Derivative variable names use prefixes z/r/t (ZRT) or z/n/e (ZNE). Set this when strain, stress or rotation will be computed later.
return_result – If true, read the generated NetCDF file with pygrt.utils.read_nc_variables().

返回:

The synthesized NetCDF data when return_result is true; otherwise None.

compute_static_syn(*args, **kwargs)[源代码]: Legacy interface renamed to static_syn(); calling it raises an error.