SpectralParameters/src/SpectralAnalysisApp.py

from obspy import read, read_events, read_inventory

from SpectralAnalysis import SpectralAnalysis
from base_logger import getDefaultLogger


def count_matching_stations(stream, inventory):
    """
    Count how many stations from the given stream are present in the inventory.

    Parameters:
    - stream: ObsPy Stream object containing seismic data.
    - inventory: ObsPy Inventory object containing network and station metadata.

    Returns:
    - count: Number of stations from the stream found in the inventory.
    """
    # Extract unique station codes from the stream
    stream_stations = set([trace.stats.station for trace in stream])

    # Extract unique station codes from the inventory
    inventory_stations = set([station.code for network in inventory for station in network.stations])

    # Count how many stations from the stream are in the inventory
    count = len(stream_stations.intersection(inventory_stations))

    return count, len(stream_stations)

def filter_stream_by_picks(st, qml):
    """
    Select from stream only traces that have an associated pick.

    Parameters:
    - st: ObsPy Stream object containing seismic data.
    - qml: ObsPy Catalog object containing phase picks.

    Returns:
    - st2: ObsPy Stream object containing seismic data.
    """
    st2 = read().clear() #create blank stream object

    for pick in qml.picks:
        net = pick.waveform_id.network_code
        sta = pick.waveform_id.station_code
        loc = pick.waveform_id.location_code
        cha = pick.waveform_id.channel_code[:2] + "*" #wildcard to obtain all 3 components
        NSLC = net + '.' + sta + '.' + loc + '.' + cha
        st2 = st2 + st.select(id=NSLC) #add trace to new stream

    st2.merge() #remove duplicate traces
    return st2


def main(waveforms, network_inventory, picks_qml, event_latitude, event_longitude, event_depth, vp, vs, density, taper_type, taper_len,
         window_len, freq_min, freq_max, min_sn_ratio=1, min_energy_ratio=1, p_wave_analysis=True, s_wave_analysis=True,
         raytracing=True, allow_station_elevations=False,
         source_model="Madariaga", sp_fit_model="FBoatwright", norm="L2", z_threshold=6,
         q_min=1, q_max=400, velocity_model=None):
    """
    Main function for application: SPECTRALPARAMETERS
    Arguments:
        waveforms: path to input file of type 'seismogram'
        Velocity model: path to input file of type 'velocity_model'
        network_inventory: path to input file of type 'inventory'
        event_qml: path to input file of type 'quakeml_seismogram_picks'
        raytracing: parameter of type 'BOOLEAN'
        save plots: parameter of type 'BOOLEAN'
    Returns:
        File(s) named 'Katalog sejsmiczny' of type 'quakeml_seismogram_picks' and format 'XML' from working directory
        :param event_qml:
    """

    logger = getDefaultLogger(__name__)

    config = {
        # filenames
        "P_WAVE_SP_FITTING_RESULTS_JSON": "P_wave_analysis_sp_fit_solutions.json",
        "S_WAVE_SP_FITTING_RESULTS_JSON": "S_wave_analysis_sp_fit_solutions.json",
        "P_WAVE_JK_RESULTS": "P_wave_analysis_JK_solutions.csv",
        "S_WAVE_JK_RESULTS": "S_wave_analysis_JK_solutions.csv",
        "P_WAVE_SP_FITTING_RESULTS": "P_wave_analysis_sp_fit_solutions.csv",
        "S_WAVE_SP_FITTING_RESULTS": "S_wave_analysis_sp_fit_solutions.csv",
        "EVENT_RESULTS": "results.csv",
        "velocity_mat_file": velocity_model,
        # app options
        "P_wave_analysis": p_wave_analysis,
        "S_wave_analysis": s_wave_analysis,
        "raytrace": raytracing,
        "allow_station_elev": allow_station_elevations,

        #event location
        "longitude": event_longitude,
        "latitude": event_latitude,
        "depth": event_depth, # in meters

        # phase hints
        "p_phase_hint": ["P", "Pg"],
        "s_phase_hint": ["S", "Sg"],

        # source parameters
        "vp": vp,
        "vs": vs,
        "density": density,

        # window parameters
        "taper_len": taper_len,
        "window_len": window_len,
        "taper_type": taper_type,
        # filter parameters
        "freq_min": freq_min,
        "freq_max": freq_max,
        # frequency signal-to-noise ratio
        "min_sn_ratio": min_sn_ratio,
        # if station signal-to-noise ratio
        "min_energy_ratio": min_energy_ratio,

        # L1 / L2
        "norm": norm,
        # Spectrum fitting model: FBrune / FBoatwright
        "sp_fit_model": sp_fit_model,
        # Source model: Madariaga / Brune
        "sp_source_model": source_model,

        # outliers detection
        "z_threshold": z_threshold,

        # q - damping limits
        "q_min": q_min,
        "q_max": q_max,

    }

    # reading files
    stream = read(waveforms)
    inv = read_inventory(network_inventory)

    qml = read_events(picks_qml).events[0]

    if not stream:
        msg = "MSEED file is empty"
        logger.error(msg)
        raise Exception(msg)
    if not inv or len(inv) == 0:
        msg = "Inventory file is empty"
        logger.error(msg)
        raise Exception(msg)
    else:
        stream = filter_stream_by_picks(stream, qml)
        if len(stream) == 0:
            msg = "No waveform data was found corresponding to the provided picks"
            logger.error(msg)
            raise Exception(msg)
        else:
            count, total = count_matching_stations(stream, inv)
            logger.info(f"{count} out of {total} stations in stream are in the inventory")
            if count == 0:
                msg = "No stations from MSEED in the Network Inventory"
                logger.error(msg)
                raise Exception(msg)
            else:
                sa = SpectralAnalysis(stream=stream, event=qml, inventory=inv, config=config, logger=logger)
                sa.run()