Note
Go to the end to download the full example code.
04.1 Magnitude sensitivity for homogeneous model: test case 4#
This is an example of computing magnitude sensitivity for a homogeneous model.
Test case 4:
Geometry: combination of station geometry and two DAS geometries (2 vertical boreholes)
Velocity model: homogeneous, Vp = 4300 m/s, Vp/Vs = 1.73, Qp = 100, Qs = 50, rho = 2300 kg/m^3
Processing: magnitude sensitivity computation for individual geometries and for combined geometry (for grid slices only).
Load all necessary packages#
Input#
Read geometry#
station_geometry_file_path = '../data/homo/test_04/TMPL_net_geometryA.txt'
das_geometry_file_path_1 = '../data/homo/test_04/TMPL_net_geometryB.txt'
das_geometry_file_path_2 = '../data/homo/test_04/TMPL_net_geometryC.txt'
station_geometry_df = pndmod.io.read_geometry(station_geometry_file_path)
das_geometry_df_1 = pndmod.io.read_geometry(das_geometry_file_path_1)
das_geometry_df_2 = pndmod.io.read_geometry(das_geometry_file_path_2)
### Show
print("\nStation geometry:")
print(station_geometry_df.head())
print(station_geometry_df.attrs)
print("\nDAS 1 geometry:")
print(das_geometry_df_1.head())
print(das_geometry_df_1.attrs)
print("\nDAS 2 geometry:")
print(das_geometry_df_2.head())
print(das_geometry_df_2.attrs)
Station geometry:
Name Y X Z NoiseLevel
0 ST1 0.0 0.0 -0.0 5.000000e-09
1 ST2 0.0 5000.0 -0.0 5.000000e-09
2 ST3 4330.0 2500.0 -0.0 5.000000e-09
{'file_path': '../data/homo/test_04/TMPL_net_geometryA.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'velocity', 'noise_unit': 'm/s', '_is_combined': False}
DAS 1 geometry:
Y X Z NoiseLevel
0 0.0 2500.0 -0.0 5.000000e-12
1 0.0 2500.0 20.0 5.000000e-12
2 0.0 2500.0 40.0 5.000000e-12
3 0.0 2500.0 60.0 5.000000e-12
4 0.0 2500.0 80.0 5.000000e-12
{'file_path': '../data/homo/test_04/TMPL_net_geometryB.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}
DAS 2 geometry:
Y X Z NoiseLevel
0 2000.0 2500.0 -0.0 5.000000e-12
1 2000.0 2500.0 20.0 5.000000e-12
2 2000.0 2500.0 40.0 5.000000e-12
3 2000.0 2500.0 60.0 5.000000e-12
4 2000.0 2500.0 80.0 5.000000e-12
{'file_path': '../data/homo/test_04/TMPL_net_geometryC.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}
Combine geometry#
combined_geometry_df = pndmod.io.combine_geometry(station_geometry_df,
das_geometry_df_1,
das_geometry_df_2,
names=['Stations',
'DAS_1',
'DAS_2'])
### Show combined
print("\nCombined geometry:")
print(combined_geometry_df)
print(combined_geometry_df.attrs)
print("\nGeometry names:")
print(combined_geometry_df.get_names())
print("\nCombined status:",
pndmod.io.is_combined_geometry(combined_geometry_df))
### Show recovered
print("\nStation geometry (recovered):")
print(combined_geometry_df.recover('Stations').head())
print("\nCombined status:",
pndmod.io.is_combined_geometry(combined_geometry_df.recover('Stations')))
print("\nDAS 1 geometry (recovered):")
print(combined_geometry_df.recover('DAS_1').head())
print("\nCombined status:",
pndmod.io.is_combined_geometry(combined_geometry_df.recover('DAS_1')))
print("\nAttributes of DAS 1 (recovered):")
print(combined_geometry_df.recover('DAS_1').attrs)
print("\nDAS 2 geometry (recovered):")
print(combined_geometry_df.recover('DAS_2').head())
print("\nCombined status:",
pndmod.io.is_combined_geometry(combined_geometry_df.recover('DAS_2')))
print("\nAttributes of DAS 2 (recovered):")
print(combined_geometry_df.recover('DAS_2').attrs)
Combined geometry:
Name Y X Z NoiseLevel
source
Stations 0 ST1 0.0 0.0 -0.0 5.000000e-09
1 ST2 0.0 5000.0 -0.0 5.000000e-09
2 ST3 4330.0 2500.0 -0.0 5.000000e-09
DAS_1 0 NaN 0.0 2500.0 -0.0 5.000000e-12
1 NaN 0.0 2500.0 20.0 5.000000e-12
... ... ... ... ... ...
DAS_2 46 NaN 2000.0 2500.0 920.0 5.000000e-12
47 NaN 2000.0 2500.0 940.0 5.000000e-12
48 NaN 2000.0 2500.0 960.0 5.000000e-12
49 NaN 2000.0 2500.0 980.0 5.000000e-12
50 NaN 2000.0 2500.0 1000.0 5.000000e-12
[105 rows x 5 columns]
{'original_attributes': {'Stations': {'file_path': '../data/homo/test_04/TMPL_net_geometryA.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'velocity', 'noise_unit': 'm/s', '_is_combined': False}, 'DAS_1': {'file_path': '../data/homo/test_04/TMPL_net_geometryB.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}, 'DAS_2': {'file_path': '../data/homo/test_04/TMPL_net_geometryC.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}}, '_is_combined': True}
Geometry names:
['Stations', 'DAS_1', 'DAS_2']
Combined status: True
Station geometry (recovered):
Name Y X Z NoiseLevel
0 ST1 0.0 0.0 -0.0 5.000000e-09
1 ST2 0.0 5000.0 -0.0 5.000000e-09
2 ST3 4330.0 2500.0 -0.0 5.000000e-09
Combined status: False
DAS 1 geometry (recovered):
Y X Z NoiseLevel
0 0.0 2500.0 -0.0 5.000000e-12
1 0.0 2500.0 20.0 5.000000e-12
2 0.0 2500.0 40.0 5.000000e-12
3 0.0 2500.0 60.0 5.000000e-12
4 0.0 2500.0 80.0 5.000000e-12
Combined status: False
Attributes of DAS 1 (recovered):
{'file_path': '../data/homo/test_04/TMPL_net_geometryB.txt', 'original_coord_units': 'm', 'coordinate_system': 'local', 'noise_type': 'strain', 'noise_unit': 'ε', 'gauge_length': 10.0, '_is_combined': False}
DAS 2 geometry (recovered):
Y X Z NoiseLevel
0 2000.0 2500.0 -0.0 5.000000e-12
1 2000.0 2500.0 20.0 5.000000e-12
2 2000.0 2500.0 40.0 5.000000e-12
3 2000.0 2500.0 60.0 5.000000e-12
4 2000.0 2500.0 80.0 5.000000e-12
Combined status: False
Attributes of DAS 2 (recovered):
{'file_path': '../data/homo/test_04/TMPL_net_geometryC.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_04/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(-3000, 8000 + 250, 250)
gy = np.arange(-3000, 8000 + 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("grid_points.shape:",grid_points.shape)
Number of grid points: 40500
grid_points.shape: (40500, 3)
Plot geometry and imaging grid#
_ = pndvis.geometry.plot_stations_view(df=combined_geometry_df,
grid_points=grid_points,
plot_names=False,
grid_marker_size = 1,
station_marker_size=100)
_ = pndvis.geometry.plot_stations_view(df=combined_geometry_df,
grid_points=grid_points,
plot_names=False,
grid_marker_size = 1,
station_marker_size=1,
view = 'xz')
_ = pndvis.geometry.plot_stations_view(df=combined_geometry_df,
grid_points=grid_points,
plot_names=False,
grid_marker_size = 1,
station_marker_size=1,
view = 'yz')
Prepare to compute magnitude sensitivity#
# Get medium parameters
density = velocity_df.at[0,'Rho']
v_p = velocity_df.at[0,'Vp']
Q_p = velocity_df.at[0,'Qp']
v_s = velocity_df.at[0,'Vp']/velocity_df.at[0,'VpVsRatio']
Q_s = velocity_df.at[0,'Qs']
print("density, v_p, Q_p, v_s, Q_s:",density,v_p,Q_p,v_s,Q_s)
# Specify source paramerers
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)
density, v_p, Q_p, v_s, Q_s: 2300 4300 100 2485.549132947977 50
Compute and plot magnitude sensitivity for independent geometries#
print("Computing and plotting magnitude sensitivity for independent geometries...")
# Loop over geometries
for geometry in combined_geometry_df.get_names():
print(f"Computing and plotting magnitude sensitivity for {geometry} geometry...")
start_time = time()
for wave_mode in ['P','S','PS']:
X_slice_mag_sens = pndmod.core.get_mag_sensitivity(grid_coords=X_plane_points,
geometry_df=combined_geometry_df.recover(geometry),
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
Y_slice_mag_sens = pndmod.core.get_mag_sensitivity(grid_coords=Y_plane_points,
geometry_df=combined_geometry_df.recover(geometry),
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
Z_slice_mag_sens = pndmod.core.get_mag_sensitivity(grid_coords=Z_plane_points,
geometry_df=combined_geometry_df.recover(geometry),
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens, yi=gy, zi=gz,
plot_title=f"Magnitude sensitivity for {geometry} geometry at X={slice_X} m: YZ view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens, xi=gx, zi=gz,
plot_title=f"Magnitude sensitivity for {geometry} geometry at Y= {slice_Y} m: XZ view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens, xi=gx, yi=gy,
plot_title=f"Magnitude sensitivity for {geometry} geometry at depth {slice_Z} m: Map view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
Computing and plotting magnitude sensitivity for independent geometries...
Computing and plotting magnitude sensitivity for Stations geometry...
Computation time: 0.21479272842407227 seconds
Computing and plotting magnitude sensitivity for DAS_1 geometry...
Computation time: 0.3876361846923828 seconds
Computing and plotting magnitude sensitivity for DAS_2 geometry...
/home/runner/work/PyNetDesign/PyNetDesign/pynetdesign/visualization/sensitivity.py:115: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
fig, ax = plt.subplots(figsize=(10, 10))
Computation time: 0.48102450370788574 seconds
Compute and plot magnitude sensitivity for combined geometry#
print("Computing and plotting magnitude sensitivity for combined geometry...")
start_time = time()
for wave_mode in ['P','S','PS']:
X_slice_mag_sens_comb = pndmod.core.get_mag_sensitivity(grid_coords=X_plane_points,
geometry_df=combined_geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
Y_slice_mag_sens_comb = pndmod.core.get_mag_sensitivity(grid_coords=Y_plane_points,
geometry_df=combined_geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
Z_slice_mag_sens_comb = pndmod.core.get_mag_sensitivity(grid_coords=Z_plane_points,
geometry_df=combined_geometry_df,
velocity_df=velocity_df,
params=params,
wave_mode=wave_mode)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=X_slice_mag_sens_comb, yi=gy, zi=gz,
plot_title=f"Magnitude sensitivity for combined geometry at X={slice_X} m: YZ view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Y_slice_mag_sens_comb, xi=gx, zi=gz,
plot_title=f"Magnitude sensitivity for combined geometry at Y= {slice_Y} m: XZ view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
_ = pndvis.sensitivity.plot_sensitivity_slice(sens_data=Z_slice_mag_sens_comb, xi=gx, yi=gy,
plot_title=f"Magnitude sensitivity for combined geometry at depth {slice_Z} m: Map view, {wave_mode}",
clb_title="$M_w$",
geometry_df=combined_geometry_df.recover(geometry),
plot_names = True if 'Name' in combined_geometry_df.recover(geometry).columns else False,
station_marker_size = 100 if 'gauge_length' not in combined_geometry_df.recover(geometry).attrs else 10)
end_time = time()
print(f"Computation time: {end_time - start_time} seconds")
Computing and plotting magnitude sensitivity for combined geometry...
Computation time: 0.7384374141693115 seconds
Total running time of the script: (0 minutes 6.533 seconds)