Note
Go to the end to download the full example code.
03.1 Magnitude sensitivity for homogeneous model: test case 3#
This is an example of computing magnitude sensitivity for a homogeneous model.
Test case 3:
Geometry: simple DAS geometry with 1 vertical borehole
Velocity model: homogeneous, Vp = 4300 m/s, Vp/Vs = 1.73, Qp = 100, Qs = 50, rho = 2300 kg/m^3
Processing: magnitude sensitivity computation with and without directionality correction (grid slices).
Load all necessary packages#
Input#
Read geometry#
geometry_file_path = '../data/homo/test_03/TMPL_net_geometryA.txt'
geometry_df = pndmod.io.read_geometry(geometry_file_path)
print(geometry_df)
print(geometry_df.attrs)
Y X Z NoiseLevel
0 0.0 0.0 -0.0 5.000000e-09
1 0.0 0.0 10.0 5.000000e-09
2 0.0 0.0 20.0 5.000000e-09
3 0.0 0.0 30.0 5.000000e-09
4 0.0 0.0 40.0 5.000000e-09
.. ... ... ... ...
96 0.0 0.0 960.0 5.000000e-09
97 0.0 0.0 970.0 5.000000e-09
98 0.0 0.0 980.0 5.000000e-09
99 0.0 0.0 990.0 5.000000e-09
100 0.0 0.0 1000.0 5.000000e-09
[101 rows x 4 columns]
{'file_path': '../data/homo/test_03/TMPL_net_geometryA.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}
Read velocity model#
velocity_file_path = '../data/homo/test_03/TMPL_velocity_model.txt'
velocity_df = pndmod.io.read_velocity(velocity_file_path)
print(velocity_df)
print(velocity_df.attrs)
Depth Vp VpVsRatio Qp Qs Rho Ep Dt Gm
0 40000 4300 1.73 100 50 2300 0.0 0.0 0.0
{'model_type': 'homo', 'model_anisotropy': 'isotropic', 'parameter_descriptions': {'Depth': 'Depth', 'Vp': 'P-wave velocity', 'Vs': 'S-wave velocity', 'VpVsRatio': 'Vp/Vs', 'Qp': 'P-wave quality factor', 'Qs': 'S-wave quality factor', 'Rho': 'density', 'Ep': 'Thomsen epsilon parameter', 'Dt': 'Thomsen delta parameter', 'Gm': 'Thomsen gamma parameter'}, 'parameter_units': {'Depth': 'm', 'Vp': 'm/s', 'Vs': 'm/s', 'VpVsRatio': '', 'Qp': '', 'Qs': '', 'Rho': 'kg/m^3', 'Ep': '', 'Dt': '', 'Gm': ''}}
Set imaging grid parameters#
# Imaging grid parameters
gx = np.arange(-7000, 7000 + 250, 250)
gy = np.arange(-7000, 7000 + 250, 250)
gz = np.arange(250, 5000 + 250, 250)
# Generate grid points
grid_points = pndmod.utils.generate_grid(x=gx,y=gy,z=gz)
print("Number of grid points:",grid_points.shape[0])
print("Number of grid points in x:",gx.size)
print("Number of grid points in y:",gy.size)
print("Number of grid points in z:",gz.size)
print("grid_points.shape:",grid_points.shape)
Number of grid points: 64980
Number of grid points in x: 57
Number of grid points in y: 57
Number of grid points in z: 20
grid_points.shape: (64980, 3)
Plot geometry and imaging grid#
_ = pndvis.geometry.plot_stations_view(geometry_df,
grid_points=grid_points,
grid_marker_size = 1,
plot_names=False,
station_marker_size=100)
_ = pndvis.geometry.plot_stations_view(geometry_df,
view = 'xz',
grid_points=grid_points,
grid_marker_size = 1,
plot_names=False,
station_marker_size=1)
_ = pndvis.geometry.plot_stations_view(geometry_df,
view = 'yz',
grid_points=grid_points,
grid_marker_size = 1,
plot_names=False,
station_marker_size=1)
Prepare to compute magnitude sensitivity#
# Specify source parameters
frequency = 10 # Hz
# Minimum S/N for detection on individual channel
min_SNR_p = 2
min_SNR_s = 2
# Minimum stations on which an event must be detected
min_stations_p = 3
min_stations_s = 3
# Save parameters
params = pndmod.classes.DetectionParameters(
f_p=frequency,
f_s=frequency,
min_SNR_p=min_SNR_p,
min_SNR_s=min_SNR_s,
min_stations_p=min_stations_p,
min_stations_s=min_stations_s,
)
# Define slices to use
slice_X = 250
slice_Y = 250
slice_Z = 2000
# Get slice planes
X_plane_points, X_plane_indices = pndmod.utils.get_slice_coordinates(grid_points=grid_points,slice_dimension=0,slice_coordinate=slice_X)
Y_plane_points, Y_plane_indices = pndmod.utils.get_slice_coordinates(grid_points=grid_points,slice_dimension=1,slice_coordinate=slice_Y)
Z_plane_points, Z_plane_indices = pndmod.utils.get_slice_coordinates(grid_points=grid_points,slice_dimension=2,slice_coordinate=slice_Z)
Retrieve minimum displacement amplitudes#
start_time = time()
min_amps_p, min_amps_s = pndmod.utils.retrieve_min_amps(geometry_df = geometry_df,
f_p = frequency,
f_s = frequency,
SNR_p = min_SNR_p,
SNR_s = min_SNR_s,
sample_interval=geometry_df.attrs['gauge_length'])
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
Computation time: 0.00017786026000976562 seconds
Compute magnitude sensitivity without directionality correction#
# Compute magnnitude sensitivity for the planes for P
print("Computing magnitude sensitivity on planes for P wave...")
start_time = time()
X_slice_mag_sens_stations_p = pndmod.homo.mag_sensitivity_grid(grid_coords=X_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_stations_p=min_stations_p,
f_p=frequency,
wave_mode='P')
Y_slice_mag_sens_stations_p = pndmod.homo.mag_sensitivity_grid(grid_coords=Y_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_stations_p=min_stations_p,
f_p=frequency,
wave_mode='P')
Z_slice_mag_sens_stations_p = pndmod.homo.mag_sensitivity_grid(grid_coords=Z_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_stations_p=min_stations_p,
f_p=frequency,
wave_mode='P')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
# Compute magnnitude sensitivity for the planes for S
print("Computing magnitude sensitivity on planes for S wave...")
start_time = time()
X_slice_mag_sens_stations_s = pndmod.homo.mag_sensitivity_grid(grid_coords=X_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_s=min_amps_s,
min_stations_s=min_stations_s,
f_s=frequency,
wave_mode='S')
Y_slice_mag_sens_stations_s = pndmod.homo.mag_sensitivity_grid(grid_coords=Y_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_s=min_amps_s,
min_stations_s=min_stations_s,
f_s=frequency,
wave_mode='S')
Z_slice_mag_sens_stations_s = pndmod.homo.mag_sensitivity_grid(grid_coords=Z_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_s=min_amps_s,
min_stations_s=min_stations_s,
f_s=frequency,
wave_mode='S')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
# Compute magnnitude sensitivity for the planes for P and S
print("Computing magnitude sensitivity on planes for P and S waves...")
start_time = time()
X_slice_mag_sens_stations_ps = pndmod.homo.mag_sensitivity_grid(grid_coords=X_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_amps_s=min_amps_s,
min_stations_p=min_stations_p,
min_stations_s=min_stations_s,
f_p=frequency,
f_s=frequency,
wave_mode='PS')
Y_slice_mag_sens_stations_ps = pndmod.homo.mag_sensitivity_grid(grid_coords=Y_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_amps_s=min_amps_s,
min_stations_p=min_stations_p,
min_stations_s=min_stations_s,
f_p=frequency,
f_s=frequency,
wave_mode='PS')
Z_slice_mag_sens_stations_ps = pndmod.homo.mag_sensitivity_grid(grid_coords=Z_plane_points,
velocity_df=velocity_df,
geometry_df=geometry_df,
min_amps_p=min_amps_p,
min_amps_s=min_amps_s,
min_stations_p=min_stations_p,
min_stations_s=min_stations_s,
f_p=frequency,
f_s=frequency,
wave_mode='PS')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
Computing magnitude sensitivity on planes for P wave...
Computation time: 0.02816009521484375 seconds
Computing magnitude sensitivity on planes for S wave...
Computation time: 0.21951770782470703 seconds
Computing magnitude sensitivity on planes for P and S waves...
Computation time: 0.24031352996826172 seconds
Plot magnitude sensitivity without directionality correction#
Plot magnitude sensitivity for P wave#
clb_range = None
plt_levels = None
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_p, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_p, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_p, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Plot magnitude sensitivity for S wave#
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_s, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_s, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_s, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Plot magnitude sensitivity for P and S waves#
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_ps, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_ps, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_ps, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Compute magnitude sensitivity with directionality correction#
# Compute magnnitude sensitivity for the planes for P
print("Computing magnitude sensitivity on planes with correction for station directionality for P wave...")
start_time = time()
X_slice_mag_sens_stations_st_dir_p = pndmod.core.get_mag_sensitivity(grid_coords=X_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='P')
Y_slice_mag_sens_stations_st_dir_p = pndmod.core.get_mag_sensitivity(grid_coords=Y_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='P')
Z_slice_mag_sens_stations_st_dir_p = pndmod.core.get_mag_sensitivity(grid_coords=Z_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='P')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
# Compute magnnitude sensitivity for the planes for S
print("Computing magnitude sensitivity on planes with correction for station directionality for S wave...")
start_time = time()
X_slice_mag_sens_stations_st_dir_s = pndmod.core.get_mag_sensitivity(grid_coords=X_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='S')
Y_slice_mag_sens_stations_st_dir_s = pndmod.core.get_mag_sensitivity(grid_coords=Y_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='S')
Z_slice_mag_sens_stations_st_dir_s = pndmod.core.get_mag_sensitivity(grid_coords=Z_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='S')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
# Compute magnnitude sensitivity for the planes for S
print("Computing magnitude sensitivity on planes with correction for station directionality for P and S waves...")
start_time = time()
X_slice_mag_sens_stations_st_dir_ps = pndmod.core.get_mag_sensitivity(grid_coords=X_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='PS')
Y_slice_mag_sens_stations_st_dir_ps = pndmod.core.get_mag_sensitivity(grid_coords=Y_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='PS')
Z_slice_mag_sens_stations_st_dir_ps = pndmod.core.get_mag_sensitivity(grid_coords=Z_plane_points,
geometry_df=geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode='PS')
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
Computing magnitude sensitivity on planes with correction for station directionality for P wave...
Computation time: 0.028563976287841797 seconds
Computing magnitude sensitivity on planes with correction for station directionality for S wave...
Computation time: 0.23161864280700684 seconds
Computing magnitude sensitivity on planes with correction for station directionality for P and S waves...
Computation time: 0.2688732147216797 seconds
Plot magnitude sensitivity with directionality correction#
Plot magnitude sensitivity for P wave#
clb_range = None
plt_levels = None
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_st_dir_p, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_st_dir_p, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_st_dir_p, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, P-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Plot magnitude sensitivity for S wave#
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_st_dir_s, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_st_dir_s, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_st_dir_s, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, S-wave",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Plot magnitude sensitivity for P and S waves#
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_stations_st_dir_ps, yi=gy, zi=gz,
plot_title="Magnitude sensitivity at X=" + str(slice_X) + " m: YZ view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_stations_st_dir_ps, xi=gx, zi=gz,
plot_title="Magnitude sensitivity at Y=" + str(slice_Y) + " m: XZ view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=1,
plot_names=False)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_stations_st_dir_ps, xi=gx, yi=gy, geometry_df=geometry_df,
plot_title="Magnitude sensitivity at depth " + str(slice_Z) + " m: Map view, P- and S-waves",
clb_title="$M_w$",
clb_range=clb_range,
plt_levels=plt_levels,
station_marker_size=10,
plot_names=False)
Total running time of the script: (0 minutes 4.116 seconds)