115 lines
4.2 KiB
Python
115 lines
4.2 KiB
Python
|
import logging
|
||
|
|
import sys
|
||
|
|
|
||
|
|
import pandas as pd
|
||
|
|
|
||
|
|
from obspy import read_inventory, Catalog
|
||
|
|
from obspy.core.inventory import Inventory
|
||
|
|
|
||
|
|
from pyoctoAssociator import PhaseAssociation
|
||
|
|
from catalogConverter import CatalogConverter
|
||
|
|
from qml2picks import qml2picks
|
||
|
|
from count_picks import count_picks_from_stations
|
||
|
|
|
||
|
|
|
||
|
|
def main(qml_picks, inventory, velocity_model, lat_min, lat_max, lon_min, lon_max, z_min,
|
||
|
|
z_max, grid_hor_extent, grid_spacing, p_vel_const, s_vel_const,
|
||
|
|
n_picks=4, n_p_picks=0, n_s_picks=0, n_p_s_picks=0, time_before=10.0, tolerance=1.0, assoc_cut_dist=250.0,
|
||
|
|
pick_match_tolerance=0.5, p_hint="P", s_hint="S"):
|
||
|
|
"""
|
||
|
|
Main function for application: PHASEASSOCIATIONAPP
|
||
|
|
Arguments:
|
||
|
|
velocity_model: input file of type 'mat'
|
||
|
|
qml_picks: input file of type 'quakeml_seismogram_picks'
|
||
|
|
inventory: input file of type 'inventory'
|
||
|
|
grid_hor_extent: parameter of type 'DOUBLE'
|
||
|
|
grid_spacing: parameter of type 'DOUBLE'
|
||
|
|
lat_min: parameter of type 'DOUBLE'
|
||
|
|
lat_max: parameter of type 'DOUBLE'
|
||
|
|
lon_min: parameter of type 'DOUBLE'
|
||
|
|
lon_max: parameter of type 'DOUBLE'
|
||
|
|
z_min: parameter of type 'DOUBLE'
|
||
|
|
z_max: parameter of type 'DOUBLE'
|
||
|
|
n_picks: parameter of type 'INTEGER'
|
||
|
|
n_p_picks: parameter of type 'INTEGER'
|
||
|
|
n_s_picks: parameter of type 'INTEGER'
|
||
|
|
n_p_s_picks: parameter of type 'INTEGER'
|
||
|
|
pick_match_tolerance: parameter of type 'DOUBLE'
|
||
|
|
p_vel_const: parameter of type 'DOUBLE'
|
||
|
|
s_vel_const: parameter of type 'DOUBLE'
|
||
|
|
p_hint: parameter of type 'STRING'
|
||
|
|
s_hint: parameter of type 'STRING'
|
||
|
|
time_before: parameter of type 'DOUBLE'
|
||
|
|
tolerance: parameter of type 'DOUBLE'
|
||
|
|
assoc_cut_dist: parameter of type 'DOUBLE'
|
||
|
|
Returns:
|
||
|
|
File(s) named 'out_catalog' of type 'data_catalog_csv' and format 'PLAIN_TEXT' from working directory
|
||
|
|
"""
|
||
|
|
|
||
|
|
inv = Inventory()
|
||
|
|
for inventory_item in inventory:
|
||
|
|
inv.extend(read_inventory(inventory_item))
|
||
|
|
|
||
|
|
if count_picks_from_stations(quakeml_file=qml_picks, inventory=inv) == 0:
|
||
|
|
logging.error(f"There are no stations in the Station Inventory file(s) corresponding to the picked times. "
|
||
|
|
f"Check if the file(s) is/are correct")
|
||
|
|
sys.exit(1)
|
||
|
|
|
||
|
|
#convert input velocity from m/s to km/s
|
||
|
|
if p_vel_const is not None:
|
||
|
|
p_vel_const = p_vel_const / 1000
|
||
|
|
if s_vel_const is not None:
|
||
|
|
s_vel_const = s_vel_const / 1000
|
||
|
|
|
||
|
|
config_pyocto = {
|
||
|
|
"velocity": {
|
||
|
|
"path": velocity_model,
|
||
|
|
"path_out": "velocity_model_1D",
|
||
|
|
"0D_P_vel": p_vel_const,
|
||
|
|
"0D_S_vel": s_vel_const,
|
||
|
|
"grid_hor_extent": grid_hor_extent,
|
||
|
|
"grid_ver_extent": z_max-z_min,
|
||
|
|
"grid_spacing": grid_spacing
|
||
|
|
},
|
||
|
|
"grid": {
|
||
|
|
"lat_min": lat_min,
|
||
|
|
"lat_max": lat_max,
|
||
|
|
"lon_min": lon_min,
|
||
|
|
"lon_max": lon_max,
|
||
|
|
"z_min": z_min,
|
||
|
|
"z_max": z_max
|
||
|
|
},
|
||
|
|
"params": {
|
||
|
|
"time_before": time_before,
|
||
|
|
"n_picks": n_picks,
|
||
|
|
"n_p_picks": n_p_picks,
|
||
|
|
"n_s_picks": n_s_picks,
|
||
|
|
"n_p_and_s_picks": n_p_s_picks,
|
||
|
|
"tolerance": tolerance,
|
||
|
|
"association_cutoff_distance": assoc_cut_dist,
|
||
|
|
"pick_match_tolerance": pick_match_tolerance
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
config_conv = {
|
||
|
|
"P_hint": p_hint,
|
||
|
|
"S_hint": s_hint,
|
||
|
|
}
|
||
|
|
picks = qml2picks(qml_picks, p_hint, s_hint)
|
||
|
|
if not picks:
|
||
|
|
# Saving empty qml file
|
||
|
|
cat = Catalog()
|
||
|
|
cat.write("detections.xml", format="QUAKEML")
|
||
|
|
else:
|
||
|
|
pa = PhaseAssociation(inventory=inv, config=config_pyocto)
|
||
|
|
events, assignments = pa.associate_phases(picks)
|
||
|
|
if not events.empty:
|
||
|
|
events['datetime'] = pd.to_datetime(events['time'], unit='s')
|
||
|
|
cc = CatalogConverter(config_conv, picks, events, assignments, "PyOcto")
|
||
|
|
cc.catalog2obspy()
|
||
|
|
cc.save_catalog_to_file("detections.xml")
|
||
|
|
else:
|
||
|
|
# Saving empty qml file
|
||
|
|
cat = Catalog()
|
||
|
|
cat.write("detections.xml", format="QUAKEML")
|